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ABSTRACT 


Numerical methods of integration of the equations of motion of a 
controlled satellite under the influence of gravity-gradient torque 
are considered. The results of computer experimentation using a number 
of Runge-Kutta, multi-step, and extrapolation methods for the numeri- 
cal integration of this differential system are presented, and par- 
ticularly efficient methods are noted. A large bibliography of 
numerical methods for initial value problems for ordinary differential 
equations is presented, and a compilation of Runge-Kutta and multi- 
step formulas is given. Less common numerical integration techniques 
from the literature are noted for further consideration. 


This report was prepared by Department of Aerophysics and Aerospace 
Engineering, Mississippi State University under Contract NAS8-28833 for 
the George C. Marshall Space Plight Center of the National Aeronautics 
and Space Administration. 



INTRODUCTION 


The integration of the equations of motion describing the dynamics 
of a body in orbit, as affected by various perturbing forces and the 
corrective action of control systems, normally requires the use of 
hybrid computers because of the difficulty and time involved in inte- 
gration of high frequency components by present numerical techniques. 

It is desirable, however, to be able to describe the vehicular 
dynamics and the vehicle response to control action numerically in order 
to take advantage of the advanced development and ease of use of digital 
computers, particularly on board the vehicle. This can only be achieved 
by improved mathematical analysis of the numerical integration of the 
initial value problem with simultaneous ordinary differential equations. 
While the numerical techniques that have become standard are adequate 
for the integration of such systems when only low frequency modes are 
involved, the efficient analysis of high frequency modes requires the 
development of new approaches. 

It is therefore necessary to develop further the mathematical anal- 
ysis of the numerical integration of systems of simultaneous ordinary 
differential equations as an initial value problem with specific ap- 
plication to the equations describing the vehicular dynamics of a con- 
trolled body in orbit. The ultimate goal of the present project is 
to develop particularly efficient numerical integration schemes for 
the equations of motion of a flexible orbiting body rotating under the 
influence of gravity gradient and control torques. 

In the present investigation an extensive bibliography of papers 
dealing with the numerical solution of systems of ordinary differential 



2 


equations has been compiled and is included in this report. Many dif- 
ferent methods of solution have been obtained therefrom to be compared. 
In the present effort comparisons of a large number of Runge-Kutta, 
multi-step, hybrid, and extrapolation methods have been made using 
the equations of motion of a rigid satellite in circular orbit ro- 
tating under the influence of gravity gradient and control torques to 
a fixed attitude, and the results of these comparisons are reported 
herein. Further effort is required to extend the comparison to other 
types of methods, and to compare <.11 methods in regard to the equations 
of motion of non-rigid satellites. 
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CHAPTER I 

EQUATIONS OF MOTION 

For a rigid space vehicle in orbit, six ordinary differential 
equations are required to specify its spatial motion — three equations 
for rotations and three for displacements. If the vehicle is con- 
sidered non-rigid, more equations (partial differential equations) are 
required to specify the deformations of the body. 

In this study the body is considered rigid. The equations of 
motion are formulated in terms of quaternions, these being functions 
of the direction cosines of the body axes. The introduction of these 
quaternions results in seven first order ordinary differential equa- 
tions which specify the motion of the body. 

The vehicle Is in orbit and the motion is assumed to be affected 
only by earth’s gravitational field. The space vehicle is to be kept 
attitude fixed, that is, its direction should be invariant for all time. 
Any deviations from this fixed direction are to be corrected by an 
on-board control system. Deviations from the fixed attitude occur 
because of the gravity gradient and imperfections in the launch 
process. These deviations are indicated by the quaternions. 


The seven first-order ordinary differential equations, in matrix 
form are: 

♦ 3GM 

I*w + w x I*w “ r x I*r + T 1.1 (3 equations) 

" "" ir| 5 " C 

3. = - _n(w)_! 


1.2 (4 equations) 



4 



is the angular velocity vector. 


w is the time rate of change of w, 


h\ 


i = 


W 


/e^sin i|)/2\ 


e 2 sin ift/2 


e^sin ^/2 
y cos ip/2 J 


is the quaternion vector, 


^ is the time rate of change of 

^ is the angle of rotation of the body’s rotation vector j 2 from a fixed 
attitude, 

is the control torque vector, a function of w and £, 



fhi 

1 12 

1 13 \ 


1 21 

±22 

1 23 


^31 

*32 

1 33j 


is the moment of inertia matrix^ 


r is the radius vector of the orbit referred to the vehicle fixed 
axes and is related to the inertial axes, fixed to the center of 
the earth, by a transformation matrix I), 

G is the universal gravitational constant, 

M is the mass of the earth, 

£2(w) is an asymmetric matrix function of w as defined below: 
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/ 


-W A 


w 2 


n(w) = 




"Wr 


-W, 


Wi 


w r 


w. 


V 


•Wj -w 


2 -w 3 


Since I is non-singular, we may write (1.1) as 


w = I ^ r x !_■*. “ w x J>w + T „J 




1.3 


Evaluating the matrix cross products and defining an asymmetric 
matrix function F(a), where a is any vector, by 


0 - a ^ <* 2 


F(a) = 


“3 0 ‘“I 


'■“2 “l 0 


we have 


_-l r 3GM 
w = I [ 


F(r) I r - F(w)I w + T ] 


1,4 


i = - jB (w) s. 


1.5 


The transformation of r^ (inertial axes) to _r (vehicle axes) is 


given by 


/cos Wgt 


—I = E o 


sin w^t 


1.6 


V 


o 
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where 

r, 2 2 2 2 

D 11 = q l ' q 2 " q 3 + q 4 
D 12 = 2< - q l q 2 + q 3 q 4' ) 

D 13 = 2 ^ q l q 3 “ q 2 q 4^ 

D 21 = 2(q l q 2 - q 3 q 4 ) 

„ 2,2 2 2 

°22 " - q l + q 2 " q 3 + q 4 

°23 = 2 ^ q 2 q 3 + q l q 4^ 

D 31 = 2(q l q 3 + q 2 q 4 } 

D 32 = 2(q 2 q 3 ” q l q 4 ) 

7 2 2 2 

D 33 * - q ! - q 2 + q 3 + q 4 

The control torque vector is given by 

T = A q + b w 
c — 

where 


1.7 


1.8 


1.9 


A is a control matrix (with dimensions of torque) , taken in the 
present comparisons as 

/4250 0 0 

A = 2 0 39,950 0 I nt - m 

\o 0 39 ’ 950 - 

b is a control matrix (with dimensions of angular momentum) : 


[B = W() b] 


b = 


-3400 0 0 

0 -8800 0 
^ 0 0 -8800 


> 


kg - m /sec 
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For the spacecraft under consideration in these comparisons, the 
moments of inertia in MKS units are: 



0 

88,950 

0 



kg-m 


2 


^ i 3 2 

with: GM - 3.98602 x 10"^ in /sec 

R 0 = 6.6525535 x 10 7 m for a 90 minute circular orbital period. 
Equations 1.4 and 1.5 can be expressed in vector form as an 
initial value problem as follows: 


dY 


1.10 


where: 


X<V - h 


Y = 


/ 


\ 




w r 


w. 



is the initial time and 



q 2 (V 
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CHAPTER II 

NUMERICAL EXPERIMENTATION PROCEDURE 


The complexities of the system preclude an analytical solution, 
and hence as the true solution is unknown, the various numerical in- 
tegration methods could only be compared with each other. The basis 
of comparison was the estimated truncation error as detailed below. 

No account of round-off error was made because this error depends on 
the computer used and the number of computations made in each method. 

We assume that the local truncation error for a p-order single- 
step method has the form 

T(t, h) = g h P+1 2.1 

where g is the principal error function (assumed essentially constant). 
Then proceeding from (t n _^, to ^n+l ■ W’ using two steps, 

each of size h, the truncation error is approximately 

T(t n+1’ h) = S 2ChP+1) 2,2 

so that 

Y «W - Y nS ■ ! * /1 2 ' 3 

With the same method, going from (t n _^, Y -^) to (t^ | in one 
step of size 2h, the truncation error is 

T(t n+1’ h) = S (2h) P "*" 1 2.4 


- 2P+1 * ^ 


Subtracting 2.3 from 2.5 gives 


2.5 


Y (h) (2h) 

n+1 n+1 


g hP +1 [2 P+1 - 2] 


2.6 
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Then approximately 


T(t n+ i, h) 


Y (h) (2h) 

n+1 n+1 

?P +1 - 2 


2.7 


Equation 2.7 was used to compare different methods. This is not 
an exact expression for the truncation error because g is not completely 
constant. Still this expression is a useful way to compare methods be- 
cause it gives some measure of the truncation error involved. This 
estimate applies strictly to single-step methods, but was used for all 
methods since it still is some measure of the error and should tend 
to zero for any method as the step size decreases to zero even though 
it is not strictly to be interpreted as truncation error for multi- 
step methods. 

Each method tested was run for the maximum step size, and then 
a number of runs were made with the step size halved successively. 
Equation 2.7 was used to estimate the error between the solution using 
step size h, and the solution using h/2. This error was calculated at 
each time step and a root mean square of these errors was calculated. 
This RMS value was used to judge different methods. 

The initial conditions used for the comparison were 



The period of orbit was 90 minutes. The maximum step size used was 
h = 0.001.* 


*A11 times and time steps without units are fractions of the orbit 
period. 
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Figure 1 illustrates a typical solution up to t = 0.25 (quarter 
of an orbit) . 

With the initial conditions used, w 2> q^, q 2 » q^ always re- 
main unchanged, while w^ and q^ oscillate. This is to be expected as 
the gravity differential on the spacecraft is the only disturbing 
moment, and hence only w^ and q^ are affected, q^ oscillates at twice 
the orbital period, and has no transient phase, w^, however, has a 
transient phase, and dictates the step size, h, necessary for stability. 
For this reason, only the truncation error committed in w^ was used to 
judge the methods tested. The error in w^ by far outweighs the error 
in q 3 . 

Most methods were unstable for h = 0.001. With each method six 
more h's were run, with h halved for each run. Figure 1 shows that 
the fast ocsillations of w^ die down at about t - 0.02. The comparison 
scheme discussed above was used to calculate the RMS truncation error 
up to t = 0.04. This interval includes a period of fast oscillation, 
up to t = 0.02, and a period of slow periodic oscillation, up to t = 
0,04, The h = 0.001 runs will obviously have the maximum truncation 
error, and each successive halving of h will reduce this error because 
as h -*■ 0 the truncation error -> 0, 

All runs were made in single precision arithmetic on a UNIVAC 


1106 computer. 
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CHAPTER III 

RESULTS OF COMPARISONS 


Runge-Kutta Methods * 

In Table I Runge-Kutta methods are compared by two parameters: 

RMS truncation error and computational time. The computational time 
is calculated as the number of function evaluations required by the 
particular method, i.e,, the number of stages, V, for Runge-Kutta 
methods. Thus computation time is in units of number of function 
evaluations. Each method compared is identified by its name and 
its equation number in Appendix I.* 

An X under a step size h, indicates that the method was unstable 
for that h. A - under a step size h, indicates that this step size 
was not run with this method. A 0 will occur for the smallest step 
size run because no truncation error can be calculated for the smallest 
step size. 

Figure 2 shows examples of the errors in w^ and q^ for various 
step sizes plotted against time t. Figure 3 shows examples of the 
RMS T plotted against step size for w^. Both these figures are for 
the fourth-order Rung e-Kutta-Rals ton method. 

Selecting a best method out of the ones tested is a judgment 
problem and depends on the users’ requirements. Most users require 
both speed and accuracy. Hence the methods of Table I were judged 
on the basis of an error-time parameter, ET . This is a judgment 


*Shahid Ahmed Siddiqi, M A Comparison of Various Order, Single- 
Step Explicit Runge-Kutta Methods Used to Solve the Equations of 
Motion of a Rigid Spacecraft in Circular Orbit M (unpublished M.S, 
thesis, Mississippi State University, 1973). 
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parameter which gives an equal weight to both the error and the compu- 
tational time of a particular method: 


ET 


, RMS 
2 P+1 


T 

2 



Based on this parameter, ET, the following observations are im- 
mediately made from Table I: 

Best fourth-order RKE: Ralston RKE(4,4)* 

Best fifth-order RKE: Butcher RKE(5,6) 

Best sixth-order RKE: Shanks RKE (6, 6) 

Best seventh-order RKE: Sarafyan RKE(7,10) 

(except for h = 0.000125) 

Best eighth-order RKE: Shanks RKE(8,10) 

(except for h = 0.000125) 

Now from among these five methods, the best choice, based on 
the ET parameter, is Shanks RKE (8, 10) with h - 0.00025. 

The results can be interpreted in another way: it is better, 

i.e., smaller ET, to use the Butcher RKE(5,6) at h = 0.000125 than 
to use the Shanks RKE(8,10) at h = 0.0005. Another observation, 
evident from Table I, is that the sixth, seventh, and eighth-order 
methods reach minimum ET's, but not at the smallest step size. 

It is again emphasized that the judgment parameter ET gives equal 
weight to error and time. 

Multi-Step Methods 

Although a complete investigation of the multi-step methods is 
desirable, the available time permitted only a limited investigation. 


*RKE(q,r) refers to an explicit Runge-Kutta method of q-order with 
r stages. 
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With the time factor in mind, three of the most promising types of 
predictors were chosen - Adams-Bashf orth (AB), Krough, and Craine- 
Klopfenstein (CK) , and three types of correctors were chosen - 
Adams-Moulton (AM), Rodabough-Wesson (RW) and Wesson. Also, a Butcher 
fifth order hybrid method was chosen to compare with the multi-step 
methods. 

The predictor and corrector equations were used to solve the 
initial value problem in various P-C combinations in both the PEC and 
PECE modes. Also the PE(CE) S mode was run with an accelerated Jacobi 
scheme; however, the acceleration parameter was found to be zero for 
this solution. Thus, the PE(CE) S mode was discarded. The result 
of the experimentation is presented in Table II. 

From Table II, the best of the three types of predictors was 
the CK predictor, and the best of the three types of correctors was 
the AM corrector. To verify these conclusions the following observa- 
tions were made: 

1) The three predictors were used in P-C combinations with the 
sixth-order AM corrector. The fourth-order CK was found 

to give a truncation error 9.4 times less than that of the 
fourth-order Krogh predictor for a time step of ,000125. 

Also, the fourth-order CK was found to give a truncation 
error 16.5 times less than that of the fourth-order AB 
predictor at the same time step. 

2) The fifth-order Butcher hybrid gave 3 times less truncation 
error than the eighth-order AB-AM combination for a time 
step of .000125. However, the sixth-order CK-AM combination 
gave 2 times less truncation error than the fifth-order 
Butcher hybrid method. 
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3) The fourth-order CK- sixth- order AM combination gave 6*5 times 
less truncation error than the eighth-order AB-AM combination, 

4) The RW correctors gave solutions that did not closely compare 
with the other Runge-Kutta and multi-step solutions. Thus, 
the RW correctors were discarded. 

5) The fourth- order CK-eighth-order AM combination had slightly 
less truncation error than the fourth order CK-eighth-order 
Wesson combination for all timesteps considered* 

Based on the merit factor, ET, introduced above, the fourth- 
order CK-eighth-order AM combination is the best of those considered. 

The best multi-step methods considered are generally less ef- 
fective than the best Runge-Kutta methods when judged by this merit 
factor, and stability limitations preclude the use of the larger 
step sizes allowed with the Runge-Kutta methods with the higher- 
order methods. It thus appears that the best Runge-Kutta methods are 
to be preferred for this system of differential equations. 

Extrapolation Methods 

Six extrapolation methods were investigated. These methods used 
rational function and polynomial function extrapolation with the 
Euler and modified midpoint algorithms. The Euler algorithm was 
used with polynomial function extrapolation with the basic time step 
being subdivided according to each of the sequences 

\ ■ V 2 “ 

and 

\ = V 2 ’ V 3 ’ V 4 ’ ••• } • 

The modified midpoint rule was used with both rational function and 
polynomial function extrapolation with each of the above sequences. 
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All methods used local extrapolation. The A-stable trapezoidal rule 
was not used since it requires global extrapolation. Further details 
on the Euler-Romberg and Bulirsch-Stoer modified midpoint method 
can be found in Reference [556]. 

Since the extrapolation methods subdivide each time step re- 
peatedly until a desired tolerance between successive extrapolations 
is obtained, comparisons with other methods is not directly possible. 
However, some conclusions can be drawn concerning the best of the 
extrapolation methods with regard to accuracy and number of function 
evaluations. 

Extrapolation methods have as variable parameters the basic step 
size and the tolerance between successive extrapolations. In this 
work there is no automatic step size correction. Results for compari- 
son purposes were obtained by making a number of computer runs at 
different basic step sizes with constant tolerance between successive 
extrapolations and then repeating the runs with successively smaller 
tolerances between successive extrapolations. These results indicate 
an optimum value of tolerance and step size for each method. As the 
basic step size increases, the number of extrapolations (function 
evaluations) needed to obtain a given tolerance decreases. However, 
for large step sizes the accuracy of the extrapolated values tends 
to decrease. The optimum values must be determined for each method 
experimentally. If the tolerance between successive extrapolations 
is too low, accuracy is poor, while if too high, instability can occur. 

All midpoint rule methods became unstable at a basic step size 
of .00037 (fraction of one orbit), some at a smaller step size. So 
apparently extrapolation cannot totally overcome the instability 
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characteristic of the midpoint rules. When using a halving sequence, 
polynomial function extrapolation required more function evaluations 
than the rational function extrapolation with little difference in 
accuracy. If a reciprocal sequence was used the polynomial function 
extrapolation required fewer function evaluations than the rational 
function extrapolation method but gave less accurate results and 
went unstable at a lower step size. Thus there is little advantage 
to either rational function or polynomial function extrapolation, and 
the reciprocal and halving sequences produced no real savings in 
function evaluations when used with the midpoint rule. 

Of the two Euler methods which both used polynomial function ex- 
trapolation, the solution obtained with the reciprocal sequence was 
very slightly more accurate than the solution obtained by successively 
halving the time step. The number of function evaluations was from 
2 to 8 times fewer for the solution with the reciprocal sequence of 
step size subdivision, depending upon the smallness of the tolerance 
between successive extrapolations. 

Partial results are presented in Figure 4 . The best two methods 
are shown, these being the Euler algorithm with polynomial function 
extrapolation utilizing a reciprocal sequence for step size subdivision 
(ERPESR) and the modified midpoint rule with rational function ex- 
trapolation utilizing a halving sequence (BSRESH) . 

The accuracy of both methods is seen to be a function of step 
size and tolerance between successive extrapolations. BSRESH ap- 
proaches the accuracy of ERPESR but requires about 40 times as many 
function evaluations at a step size of .000185 (fraction of one orbit). 
At a step size of .00037 the order of the function evaluations is the 
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same but BSRESH goes unstable at higher values of time* Extrapolation 

to higher tolerances increases accuracy , particularly at lower step 

sizes, at the expense of increased function evaluations for both 

BSRESH and ERPESH. However, there is a limit to the tolerances for 

both methods below which the solution becomes unstable. And, of 

course, if the tolerance is too large accuracy decreases, 

ERPESH had constant accuracy over a large range of tolerances 
-4 -10 

(10 - 10 ) for step sizes between ,0000231 and .000185. However, 

at larger step sizes the accuracy of the solution decreased at lower 
tolerances (10 ). 

Thus the best of the extrapolation methods is ERPESR, Euler 

algorithm with polynomial function extrapolation using a reciprocal 

sequence to subdivide the basic step size until a tolerance between 

-4 -10 

successive extrapolation of 10 to 10 is achieved. This method 
gives better accuracy than the Runge-Kutta and predictor corrector 
methods and takes about the same number of function evaluations, A 
value of angular velocity obtained by the best of the predictor 
corrector methods is shown on Figure 5 at a step size of .000185, 

This value is 1.4% lower than the value obtained by ERPESR and re- 
quired 40 function evaluations as compared to 41 function evaluations 

-4 

for ERPESR when the extrapolation was carried to a tolerance of 10 , 

or 71 function evaluation for a tolerance of 10 Thus ERPESR 

appears to require the same order of function evaluation as Runge- 
Kutta and predictor corrector methods at a better accuracy. 



19 


CHAPTER IV 

SURVEY OF TECHNIQUES FOR FURTHER CONSIDERATION 
The numerical solution of the initial value problem 

= 1 ^* 1 ) » z(o) - Z • ^.i 

has been approached in various ways as outlined below. As is noted, 
the various classes cited are not exclusive, but overlap one another. 
Quadrature Methods 

One may write formally 

t 

£(t ) = ^(t ) + / r 1 (t, y(t)) dt 4.2 

q t 

q 

In the quadrature methods, the integral in 4.2 is approximated by a 
numerical quadrature expression. For instance, suppose it is desired 
to calculate the numerical approximation, y^, at a series of times, t^, 
n=0, 1, 2,«**> called the nodes hereafter for convenience of notation, 
these nodes not necessarily being equally spaced. Then we may approxi- 
mate the integral of 4.2 by a numerical quadrature consisting of a 
linear combination of nodal values and possibly also some intermediate 
values at points interspersed among the nodes. The coefficients in the 
linear combination and the locations of the nodes and/or the inter-nodal 
points are determined by the particular type of numerical quadrature 
employed. These factors are obtained directly from quadrature expres- 
sions for some types of methods, but by only indirect reference to 
numerical quadrature in others, the former type being that strictly 
referred to as "quadrature methods 11 in many works. 
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If only previously calculated nodal values, and possibly the 
value at the node presently under consideration, are included in the 
quadrature formula, then no additional equations are required, and 
the value at the node presently under consideration is determined 
either explicitly or implicitly through the solution of a nonlinear 
equation, depending on whether or not the value at the node under 
consideration is included in the quadrature expression. Methods of 
this type are also a sub-class of the linear multi-step methods dis- 
cussed below. 

Nesterchuk {378] gives such an implicit method involving 
all previous nodal points. 

If nodal points beyond the node presently under consideration are 
included in the quadrature formula then an equation must be written for 
each of these points, and a system of equations must then be solved 
simultaneously for all the unknown nodal values involved „ This type 
method is also a sub-class of the block linear multi-step methods dis- 
cussed below. 

Day [116] gives a method of this type using Lobatto 
quadrature and successively higher moments of the differential 
equation to supply the needed equations for all the unknown 
nodal values included in the block. For linear equations 
this method achieves n+1 order with n function evaluations 
and is reported to be superior in accuracy to some Runge- 
Kutta methods. 

If inter-nodal points are involved, additional equations for 
the values at these points must be included. If these inter-nodal 
values are determined in turn by quadrature formulas, we have the 
standard Runge-Kutta methods. Within this class we have explicit 
methods if the expressions for the inter-nodal values involve only 
values at previous inter-nodal points, semi-implicit methods if these 
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expressions involve the presently considered inter-nodal value as 
well but none beyond, and implicit methods if inter-nodal values beyond 
that presently considered are included. The latter types require, 
respectively, the solution of one nonlinear equation at each inter- 
nodal point or the simultaneous solution of a system of equations, 
one for each internodal point. 

The method of Hulme [253] is a generalization of this im- 
plicit Runge-Kutta form constructed using Gaussian quadrature 
to approximate the integrals involved in a Galerkin approxi- 
mation with the solution represented by a polynomial within 
each step. This produces a piecewise continuous solution 
with essentially the same effort required by conventional 
implicit Runge-Kutta methods and the same nodal values as 
produced thereby, Axelsson [13], in a general discussion of 
quadrature methods, presents implicit Runge-Kutta methods 
based on Radau and Lobatto quadratures. Implicit Runge- 
Kutta methods are also discussed in Gourlay [194], where it 
is noted that certain methods stable for linear equations, 
i.e., with constant Jacobian matrices, may not be stable 
when the Jacobian varies as is the case with nonlinear equa- 
tions, The common trapezoidal method is a case in point, 
and a stable modification of the same order is given therefor. 

Explicit Runge-Kutta methods with provision for intrinsic 
truncation error estimation are given in Warten [536] » 

Sarafyan [448], and Zonneveld [553], Haines [220] gives a 
semi-implicit Runge-Kutta method with the coefficients chosen 
to increase stability and with the Jacobian matrix, required 
in such methods, evaluated by finite differences. Semi- 
implicit Runge-Kutta methods are also given by Allen [3], 

Sarafyan [452] fits a Hermite polynomial to the inter- 
nodal values of an explicit embedded Runge-Kutta method to 
produce a continuous approximation of order only one less 
than that of the discrete approximation. Merson [356] gives 
Runge-Kutta methods for which each inter-nodal value is 
required to be a good approximation of the solution at the 
corresponding point. Treanor [519] uses local linearization 
before the quadrature to develop a modification of the 
Runge-Kutta type. The algebraically difficult derivation 
of the coefficients involved in Runge-Kutta methods is rel- 
egated to the computer via a program of Sarafyan and Brown 
[447]* Many other Runge-Kutta methods are discussed in 
Appendix I, 
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The more general inclusion of inter-nodal points along with nodal 
points in the quadrature expression for the integral of 4.2 results 
in a combination of the above types, variously referred to as multi- 
step Runge-Kutta methods or hybrid methods, 

Rosen [430] gives an explicit multi-step Runge-Kutta 
method using at each node all the inter-nodal values used 
at the previous node, with a resultant decrease in the number 
of evaluations per step required for a given order. 

Finally, in the manner already discussed above, the inclusion of 
nodal and/or inter-nodal points beyond the node presently under con- 
sideration produces block methods of the Runge-Kutta or hybrid type, 
requiring simultaneous solution for all the unknown nodal and/or 
inter-nodal values involved. 

Block Runge-Kutta methods in which the iteration at each 
node is intentionally not carried to convergence are discussed 
by Rosser [433]. These methods require fewer evaluations per 
step for a given order than the usual point Runge-Kutta methods. 

Multi-value Methods 

At each nodal point, one may write, more generally, the value of 
the dependent variable, any of its derivatives, and any functions or 
combinations thereof as linear combinations of these quantities at 
each node and/or at inter-nodal points, and so produce the almost 
all-inclusive class of multi-value methods. 

Methods of this type are discussed in general in Gear [177]. 
Methods with variable step size are discussed in regard to 
stability in Tu [520]. 

The entire class of quadrature methods discussed above is a sub- 
class of the multi-value methods for which the value of the dependent 
variable, at only one previous nodal point is involved, and the 
linear combination includes only the first derivative, = f_ (t, %) . 
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The use of values of only the dependent variable and its first de- 
rivative in the linear combination produces the linear multi-step methods „ 
These methods are explicit if only values of the derivative at previous 
nodes are included, semi-implicit (commonly called implicit) if the de- 
rivative at the presently-considered node is included as well but none 
beyond, and implicit (commonly called "block") if nodal values beyond 
that presently under consideration are included. The latter types re- 
quire, respectively, the solution of a nonlinear equation at each node 
or the simultaneous solution of a system of equations for all the unknown 
nodal values involved. Again inter-nodal points may be included as well, 
and the so-called hybrid methods developed as discussed above. The more 
common types of predictor-corrector methods, in which values are "pre- 
dicted" by explicit multi-value forms and then "corrected" using implicit 
multi-value forms involving the predicted value as well as previously 
calculated values are in this class. Additional applications of the 
corrector may follow, using the most recently corrected values, 

Hull and Creeraer [244], in a comparison of various predictor- 
corrector forms, state that error, stability, and the order re- 
quired for a given error at least computation are less sensitive 
to order with two corrector applications than with one, with no 
significant improvement thereafter. Lambert [308] shows the equiv- 
alence of any predictor-corrector method using a finite number of 
corrector applications to some explicit multi-step method. This 
limits the expectations for stability of predictor-corrector 
methods. Donelson and Hansen [129] use a set of correctors ap- 
plied cyclically over a set of time steps and thus achieve higher- 
order stable methods than are possible with the same corrector 
used at all steps. Order 2k-l is thereby achieved for stable k- 
step methods, as opposed to the maximum order, k-1, possible with 
stability for methods using only a fixed corrector. The use of 
variable step size in predictor-corrector methods has been con- 
sidered by Van Wyk [332], and in a single-step method of the 
multi-step form by Richards, Lanning, and Torrey [419]. Other 
predictor-corrector methods are given in Appendix II. 

Many multi-step methods have been developed in which some 
coefficients are chosen to increase stability as in Krough [291], 
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with a comparison of stability plots for various methods; 

Crane and Lambert [96]; Nigro [383]; Rahme [406], with an at- 
tempt to minimize -error while maintaining stability; and Lomax 
[332], considering Runge-Kutta and hybrid methods as well* A 
multi-step method with the coefficients chosen to fit the 
characteristic roots to some of the characteristic exponentials 
of'the differential system is given by Miranker [366] for non- 
linear systems with local linearization, the method being 
implicit, and for linear systems by Liniger and Willoughby 
[323], an explicit method. Osborne [398] gives a single-step 
method with the nodes chosen to give desired characteristic 
roots. Timlake [518] increases the stability of multi-step 
methods by averaging over previous steps, determining the 
weights from the Jacobian eigenvalues. The stability of multi- 
step methods has been analyzed particularly in Karim [258,256], 
sufficient conditions for instability being given in the 
former and the effect of the predictor as well being in- 
cluded in the latter; Dahlquist [100-106], a series of classic 
papers; and Hafner [217], giving stability charts for a 
large number of Runge-Kutta and hybrid methods as well. 

Implicit multi-step methods particularly for stiff 
systems have been given \by Ratliff [411], giving a com- 
parison of several methods; Jain and Srivastava [246,247]; 

Dill [126], a systematic search for such methods; Gelinas 
[182]; Dill and Gear [125], again a search; and Gear [172]. 

Such methods achieve a particular shape of stability 
boundary suitable for stiff systems by always including 
the derivative at the present node with perhaps the deriv- 
atives at a few previous nodes. 

Tyson [521] gives an implicit multi-step method for 
nonlinear systems using local linearization about a predicted 
value, rather than the previous nodal value, that achieves 
order two greater than that of the predictor used. (Local 
linearization is used in many methods to render nonlinear 
systems tractable by methods restricted to linear systems, 
but such linearization about the previous nodal value re- 
stricts the order to two regardless of the method used.) 

Boggs [30] used Broyden iteration (involving direct approxi- 
mation of the Jacobian inverse rather than inversion of the 
Jacobian) with implicit multi-step methods. Block methods 
are given by Daniel [109,108] and by Shampine and Watts 
[469]. Methods for determining the necessary set of starting 
values (a problem equivalent to the development of block 
multi-step methods) have been given in Reimer [412], 

Rakitskii [407], and Alonzo [4]. 

The addition of inter-nodal points (hybrid methods) has 
been considered by Papian and Ball [399]; Lomax [332], who 
gives a method for choosing the coefficients in explicit 
methods of this type to improve stability; Gragg and S tetter 
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[195]; discussing the coefficients therein for maximum 
order; Gear [175]; and Hafner [217], who gives a large 
number of stability charts. Additional consideration of 
multi-step methods has been given by Spijker [492], con- 
sidering the question of direct application to higher- 
order equations or a split to a system of first-order 
equations, and by Zverskina [555], using derivatives only 
at previous nodes separated by some distance from the 
present node. 

Various other quantities, such as higher derivatives and divided 
differences, have been included in the linear combination, and methods 
using one type of quantities may be derived by transformation of 
methods using another type as in Osborne [395] and Kohfeld and Thompson 
[283], Again the above-mentioned ramifications involving inter-nodal 
points and/or nodal points beyond that presently under consideration 
may be developed. A further modification is that values at previously 
passed nodes may be subsequently changed. It is possible to derive 
methods of this type that use higher derivatives or higher-order 
divided differences, all of which are produced at each node from linear 
combinations of the same, in lieu of values at previous nodes. In 
this manner more efficient step size change can be accomplished. 

In Grobner, Kohnert, Reitberger and Wanner [207] higher 
derivatives are included in generalizations of both implicit 
Runge-Kutta and multi-step type methods. Multi-value methods 
involving higher-order derivatives are discussed in Gear 
[177,178]; Sloan [484]; Lewis and Stovall [321]; and Nordsieck 
[387], a single step form with ease of step-size change; and 
Osborne [395]. Kohfeld and Thompson [283] give a method 
using higher derivatives and inter-nodal points. The single- 
step implicit method of Urabe [522] uses second derivatives. 
Makinson [348] gives a single-step implicit method using 
higher derivatives, the form being specially constructed to 
preserve the original sparseness in the matrix solution of 
the implicit equations. Lomax [334] presents a single-step 
implicit method using second derivatives with coefficients 
chosen to increase stability. Ehle [132] and Davison [110] 
also give single-step methods using higher derivatives, and 
Reiraer [414] analyses multi-step methods using higher deriva- 
tives. 



26 


Non-Polynomial Interpolates 

In the great majority of methods the linear combinations involved 
in the multi-value methods, and thus in the quadrature methods , are 
required to be exact when the dependent variable is a polynomial of 
some degree. This restriction may be relaxed, and an even broader 
spectrum of methods may be developed by requiring these linear combina- 
tions to be exact for some other function, in particular for perhaps 
some class of functions especially adapted to the particular differential 
equation being considered. 

Loscalzo and Talbot [339] use a spline function with the 
higher derivatives matched at the nodes, which, while pro- 
ducing the same nodal values developed by a multi-step 
method, is a single-step method and gives a piecewise con- 
tinuous approximation. It does, however, require differentia- 
tion. Byrne and Chi [45] use a spline approximation of the 
intergrand in a quadrature type method requiring no differentia- 
tian. Callender [71] uses a spline approximation spanning 
several nodes to develop a block multi-step method. 

Blue and Gummel [27] and Lambert and Shaw [306] use 
rational function approximation of matrix exponentials, the 
former using coefficients chosen to increase stability and 
producing thereby Pade approximates of the exponential. 

The latter method matches higher derivatives and thus re- 
quires differentiation. Roe [84] gives a multi-step method 
using a combination of a polynomial and an exponential as 
the interpolate. Other non- polynomial interpolates have been 
used by Shaw [475]? Lambert and Shaw [303], and Pope [403], 
the first two being primarily of use when singularies are 
involved in the solution. 

Only a relatively small number of methods of this type have 
been considered, since the intention of most developments has been to 
develop methods that may be applied at least fairly well to all types 
of differential equations, hence the use of polynomial interpolation. 

In this connection the remarks of Lomax [334] are pertinent: 

"It is unlikely that ’new* combinations of linear 
equations that connect a function, u, and its derivative, 
u 1 , at a series of reference points, equispaced or not, 
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will improve existing methods for the numerical inte- 
gration of general sets of coupled ordinary differential 
equations . . ." (p. 2). M What appears to be needed are 
studies of special methods designed for special methods 
designed for special classes of equations . .. M (p, 6). 

Iterative Methods 

Methods of this class start with some approximate solution, 
such as the solution of a related but simpler differential equation or 
the result of any numerical method, and improve the approximation by 
an iterative procedure. One such method is the Lie series method. 

The extrapolation methods by which solutions with successively smaller 
truncation errors are obtained by extrapolating results with succes- 
sively smaller step-sizes may also be considered iterative in the sense 
stated. 


In reference [207] Runge-Kutta methods having error ex- 
pansions in even powers of the step size are developed to serve 
as the basis for these extrapolation methods. Lie Series 
methods have been given by Knapp and Wanner [279,280]. 

Transformation Methods 

In methods of this type the dependent and/or independent variable 

is transformed before numerical solution in order to obtain a form with 

better numerical properties. In this manner stiff systems may be 

transformed into systems with smaller ranges of eigenvalues. 

Lawson [313] gives a transformation of the dependent 
variable which reduces the stiffness of the system, but the 
method requires the use of matrix exponentials. A similar 
transformation of the dependent variable is used by Jain 
[259]. Other transformations of the dependent variable are 
given by Decell, Guseraan and Lea [117], requiring evaluation 
of the Jacobian inverse, and by Calahan [67]. 

Higher-Order Equations 

As is well known, any higher-order differential equation can be 


broken into a set of simultaneous first-order equations to which the 
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preceding methods are applicable* Generalizations of the classes cited 
follow readily, however, so that application may be made directly to 
the higher-order equation if such is thought to be preferable. This 
is still something of an open question. 

Osborne [397] gives a quadrature method for linear higher- 
order equations* Day [112,113] has implicit Runge-Kutta type 
methods using Gaussian and Lobatto quadrature, respectively, 
for second-order equations with no first derivative. Cooper 
[91,93] gives a generalization of the implicit Runge-Kutta 
methods to equations of any order and states that the split 
into a set of first-order equations decreases the local accuracy 
and increases the computation required, with little effect 
on the global accuracy. 

The multi-value methods using higher derivatives are di- 
rectly applicable to higher-order equations, and Gear [177] 
states also that direct solution of the higher-order equations 
may be faster than that of the split system of first-order 
equations. Methods of this type are also given in Gear [178, 

484 ] and Allen [2], Adrianova [1] found less error without 
splitting to the first-order system. However, Spijker [503] 
states that round-off error is reduced by splitting into first- 
order equations. 

Conclusions and Directions 

Implicit methods in general are much more stable than explicit 
methods, but may still be less efficient because of the difficulty 
of the solution of the system of simultaneous equations that is 
required at each time step. The key to the use of the more stable im- 
plicit schemes is thus the iterative solution of this system, and ef- 
fort should be directed toward the development of faster and more ap- 
propriate iterative schemes for use in these methods. There are, of 
course, a great many general iterative schemes that have never been 
applied in implicit numerical solutions of ordinary differential 
equations. 

Block methods with a set number of iterations, such as that of 
Rosser [432], can achieve a given order with fewer evaluations per step 
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than for point methods and are thus attractive. The number of evalua- 
tions required may also be decreased by the retention of previously 
calculated inter-nodal values as in the method of Rosen [430]. 

The many possible forms of the multi-value methods of Gear [177], 
in particular those that update previously calculated values, may offer 
more efficient methods than the conventional multi-step forms. Stable 
methods of higher order than possible for multi-step methods are at- 
tainable in this class. Methods of this type may be developed from the 
existing multi-step methods by linear transformations, and yet may be 
more efficient in computation. 

The local linearization of nonlinear systems about a predicted 
value, Tyson [521], rather than the value at the previous step, allows 
linear methods to be applied without limiting the resultant overall 
order for the nonlinear system to two. This procedure has received 
little attention but could widen considerably the scope of application 
of many methods applicable only to linear systems. 

The use of a cyclic set of correctors, Donelson and Hansen [129], 
rather than a fixed corrector, over a series of time steps yields 
stable methods of higher order than possible with a fixed corrector 
and this should be pursued. 

Methods using non-polynomial interpolates specialized for particu- 
lar systems of differential equations should be considered. Here speed 
may quite possibly be gained from the loss of generality. In the same 
respect, preliminary transformations of the dependent variable adapted 
to the particular differential system should also be considered. 



30 


TABLE I 


Explicit Runge-Kutta 
RKE(p.V) 

EQ. 

IN 

AP- 

PEN- 

DIX 

-1 

V 

T 

I 

M 

E 

RMS TE x 10 ~ 10 

ET - (RMS TE/2 P+ 1 -2)(V/h) x 10* 

h - 

.001 

.0005 

.00025 

.000125 

.0000625 

.00003125 

.000015625 
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A. 26 

4 

X 

1.578xl0 4 

1 . 666 * 10 3 

1.05*10 2 
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2.683 

0 



4.206-4 

8.99-5 

1.10-5 

1.349-6 

1.145-6 
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4 

X 

1.578xl0 4 

1 . 6 B 6 * 10 3 
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H 
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1.527-6 

i 
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4 

X 
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1 . 04* 10 2 

6,325 





4.173-4 

8.43-5 i 

1.11-5 

1.349-6 
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6 

X 


IP^WI 


0 


_ 



1.95x10-4 
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s 






5 



5.34X10 1 


0 






6.013-5 

2.07-6 

1.39-7 






6 

55IE5H 
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IP 
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_ 

_ 



2.876 

Htt 

4.58-6 

2, 77-7" ” 
. 




Fehlberg(5,6> 

A. 49 

6 

9 


pmtnia 


0 

_ 




8,98-5 

4.0-6 





Luther(5,6) 

A. 35 

6 

—I 



■5SSBI 

0 



Newton-Cotes Quadrature 


warn 

7.05-4 

3.50-5 

2.08-6 




Luther <5, 6 ) 

A. 43 

6 


ga^gwi 

6.90xl0 2 

2.06X10 1 

0 



Gauss Quadrature 


■ 

5.32-4 


1.59-6 




Luther (5,6) 

A. 44 

6 

■ 



8.05 

0 

IfMl 


Radau Quadrature 


■ 

1.976-4 

1.05-5 





Luther (5,6) 

A. 45 

6 


IIMJI 

2.70xl0 2 


0 

- 


Lobatto Quadrature 



1.976-4 

1.05-5 




i 

Iff 

A. 47 

5 

9.898*10 U 

1.057xl0 3 

2.72xl0 2 

1.297X10 1 

0 

mam 

■ 

SyHH wKKKKm 


7.982+3 

1.64-4 

8.78-6 

8.37-7 





A. 47 

5 

9.692X10 11 


l.lSlxlO 3 


0 


pjjpiip 




2.22-4 

3.71-5 

7.27-6 




Modified Shanks<5,4) 


L 


FBBWB 

4.20xl0 3 
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mm/m 


N-100 

HN 


KKSyHI 


2.65-5 



— 

Butcher ( 6 , 7) 


■ 

K 
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0 

Newton-Cotes Quadrature 

■ 


1.16-4 

2.63-6 

7.947-8 
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3.18-7 

■■Hi 
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6 

®SB58 
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.894 

.6325 
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4.04-7 


Sarafyan(7,10) 

A . 68 

10 


3.64+1 
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0 



Newton-Cotes Quadrature 
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The signed, number following each entry indicates multiplication by the corresponding power of 10. 
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Figure 2. Errors in and q 3 versus time, up to 1 = . 
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MULTIPLE OF BASIC TIME STEP 
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Figure 3, KMS error In versus step size. 
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Figure 4* Comparison of ERPESR and BSRESH Extrapolation Methods - Time Step Effects 
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APPENDIX I 
RUNGE-KUTTA METHODS 

This appendix is intended as a reference for various Runge-Kutta 
methods available in current literature. 

Short form notation will be used henceforth: an Explicit p order, 

V stage Runge-Kutta method is written as RKE(p,V). 

RKSI(p,V) - semi-implicit 
RKI (p , V) - implicit 
RKMS(p,V) - multi-step 

The Runge-Kutta coefficients will be presented in the following 

form: 



M. < V 
1 


In the following sections various references, where different 
Runge-Kutta methods may be found, are listed. 

Many authors have developed the equations leading to the deriva- 
tion of the coefficients defining a particular order Runge-Kutta 
method. These equations will be henceforth called the deriving equa- 
tions; and these authors will be listed first. If a reader desires 
to develop his own Runge-Kutta method, he can do so by solving the 
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deriving equations. The field is still wide open, many Runge-Kutta 
methods can still be developed. 

A. 1 Explicit Runge-Kutta Methods 
RKE(p.V) : 

Restating the equations for a RKE(p,V) 


i-1 



The deriving equations for a RKE(p,V), based on a Taylor series 
analysis, (as developed in Ch. 2) are listed in; Butcher [3], up 
to a 8th order method; and Fehlberg [4, 12], up to a 8th order method. 

A formulation for solving the non-linear deriving equations, for 
a RKE, on a computer is being developed and the principle is listed 
in Saraf yan and Brown* [ 7 ] . 

The tedious Taylor series analysis for the deriving equations 
can be replaced by a, practically, equally tedious 'Quadrative 
Method. 1 The Quadrative Method for the deriving equations is, 
however, a more convenient form for use, as it is general for any 
older method, unlike the Taylor series form. Quadrative form de- 
riving equations for any order RKE's are listed in Rosen [5]. 

An excellent reference for the general Runge-Kutta class of 
methods, all in one handy cover, is Lapidus and Seinfeld [6]. This 
book also lists commonly used Runge-Kutta methods of various orders. 
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References for RKE of specific orders are now listed. References 
for deriving equations will be given first, these will be the deriving 
equations of that specific order method. 

A. 1-1 RKE (1,1) : 

The famous (infamous?) Euler method may be considered to be such 
a method 

Euler 
A. 3 



The deriving equations, specifically for a RKE (2, 2) are listed 
in Ralston [8] and Eehlberg [12]. In Ralston [8] the free parameter 
was manipulated to give a Minimum Truncation Error Bound, henceforth 
referred to as TEB. It must be noted here that, the truncation error 
is not minimized but rather its bound is minimized. This may or may 
not result in a minimized truncation error. However, as the exact 
value of the truncation error depends on the differential equation 
being integrated, this course is the only one open for optimizing 
a method. Bounds of these kinds are quite conservative. 

Lapidus and Seinfeld [6] list 3 such methods. One of them, the 
Heun form, is Ralston's optimum. 

Heun 

(optimum) 
A. 4 
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0 

1/2 1/2 


0 1 


Mid-point Rule 
A. 5 



Euler-Cauchy 
A. 6 


^ en dx independent of Y, the Heun form becomes a 3rd order 
method. Johnston [9], Kuntzmann[10] , and King [11] also find opti- 
mum methods, these are the Heun form, eq. A. 4. 

Fehlberg [12] has developed methods of order p, coupled with a 
p+1 order method. Suitable choice of free parameters is made so as 
to minimize the leading term of the truncation error. This results 
in a larger permissable step size. Comparison of the solution of the 
p+1 order with the p order is used to control step size. 

The p+1 order shares most of the coefficients of the p order; 
hence with only 1 or 2 additional stages, a p+1 order solution is 
available. Though the Fehlberg methods require more stages than 
conventionally used; this extra calculation results in stepsize con- 
trol, which is well worth the cost. 

The automatic step size control feature of these Fehlberg 
methods make them computationally more efficient (speed and overall 
number of function evaluations), than conventional methods. This 
is so; because when step size controls, like the Richardson Extrapo- 
lation procedure (see Lapidus and Seinfeld [6]), are applied to 
conventional methods; the computational effort doubles. In the 
opinion of the author; these methods when used without the coupled 
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]j>+i order method, are still useful because the TEB has been carefully 
minimized. The p+1 order method can also be used by itself. 

Henceforth Fehlberg methods will be written as: RKE[(p,V); 


(p+1, V+n)] n = 

= 1,2,3 

• 



Fehlberg [12] RKE[(2, 

2); (3,)]: 



0 





13 

II 



1 



Fehlberg 

P=3 1/2 

1/4 1/4 



A. 7 


1/2 1/2 

- w' 

s for p=2 



1/6 1/6 

2/3 - w' 

s for p=3 


0 





1/4 

1/4 



Fehlberg 

p=2 27/40 

-189/800 

729/800 


A. 8 

p=3 1 

214/891 

1/33 

650/891 



214/891 

1/33 

650/891 

p=2 


533/2106 

0 

800/1053 -1/78 

P=3 

A. 1-3 RKE(3, 3) : 






The deriving equations for such methods are listed in Fehlberg 
[12], and in Ralston [8], as a two parameter family. In Ralston [8] 
the TEB was minimized by suitable choices of the two free parameters, 
and the 'optimum' method is: 


Ralston Optimum 
A. 9 
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Johnston [9], Kuntzm^nn [10] and King [11] also list such opti- 
mum methods. 

The two optimum methods of King [11] are, however, of a fourth 
and fifth order when is independent of Y. When ^ is dependent on 
Y his methods have a slightly larger TEB than Ralston's method (the 
method which becomes fourth order has a slightly smaller TEB than 
the method which becomes fifth order.) 



.3550510257 

.8449489743 


.3550510257 

-.4021612205 1.247110195 


.1111111111 


.5124858262 


King optimum (when 
independent of 
Y — -fourth order) 

A. 10 


.3764030627 

King optimum (when 
dY 

-jjr independent of 
Y — fifth order) 

A. 11 


0 


(10-2/13) 

6 

(10-2/13) 

(1+/13) 

6 

-..0581020 .8256939 


.2071768 .3585646 .4342585. 


Kunt zmann 
optimum 
A. 12 
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Lapidus and Seinfeld [6] list three methods: 


0 

1/2 1/2 

1 -1 

2 


Classic 
A. 13 

1/6 

2/3 

1/6 


0 



Heun 

1/3 1/3 



A. 14 

2/3 0 

2/3 



1/4 

0 

3/4 


0 




2/3 2/3 



Nystrom 

2/3 0 

2/3 


A. 15 

1/4 

3/8 

3/8 



Fehlberg [12] gives methods which incorporate automatic step 
size adjustment . 

Fehlberg [12] RKE[<3,4); (4,5)]: 
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0 




Fehlberg 

in 

2/7 



A. 17 

7/14 

77/900 

343/900 



P=3 35/38 

p=4 1 

805/1444 

79/490 

-77175/54872 

0 

97125/54872 

2175/3626 

2166/9065 


79/490 

0 

2175/3626 

2166/9065 p«3 


229/1470 

0 

1125/1813 

13718/81585 1/18 p=4 


The Fehlberg coupled second order methods, eqa. A. 7 and A. 8, 
could be used as third order methods. 


The author developed optimum methods of a special kind, these 
ethods require only two, instead of three function evaluations per step 


Back Step Method ; 


0 


-24/100 

-24/100 

76/100 

19627/10200 -475/408 


409/684 -7/36 36/57 


Butcher [18] developed such fifth order processes, using negative 
c's. This back step method has a built-in advantage; in every step 
except the first; k 3 from (T a§ 1/n) may be used as the k 2 for (T Y ) , 
This results in one less function evaluation per step, hence has V = 2 
rather than V - 3, a considerable advantage. A disadvantage of back step 
methods is that a starting method is required. For a starting method 
solve the deriving equations with c 3 o and any cholce of 
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I_I. Front Step Methods : 



These front step methods have a built in advantage: in every 

step except the last; k- from (T ' Y ) may be used as the k„ for 

n n 2 

^ T n+l’ ^n+1^ * ^ike t ^ ie hack step method, this results in a V = 2 
rather than V = 3 method. The front step method is self starting 
and so has an advantage over the back step method. 

Each of the methods eqs. A. 18, A. 19 and A. 20 was optimized using 
the procedure found in Ralston [8]. The c^s and c 3 's were inter- 
changed: a TEB was calculated using a particular c 2 and c 3 ; then an- 

other TEB was calculated using the same c 2 as a c 3 , and the same c 3 
as a c 2# Eqs. A. 18, A. 19 and A. 20 are the resulting optimum methods, 
Method (p, V) Error Bound (TEB) 


Ralston A. 9 (3,3) 

.1111 ML 3 h 4 

King A. 10 (3,3) 

.1389 ML 3 h 4 

King A. 11 (3,3) 

.1391 ML 3 h 4 

Classis A. 13 (3,3) 

.5 ML 3 h 4 

Heun A. 14 (3,3) 

.6667 ML 3 h 4 

Nystrom A. 15 (3,3) 

.25 ML 3 h 4 

Back step A. 18 (3,2) 

.3493 ML 3 h 4 

Front step A. 19 (3,2) 

.3933 ML 3 h 4 

Front step A. 19 (3,2) 

.3649 ML 3 h 4 
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A. 1-4 RKE(4, 4) ; 

More has been written and analyzed about the RKE(4,4) than any 
other Runge-Kutta method. The references listed here try to sample 
all the aspects of these methods. 

The deriving equations are listed in Butcher [3], Rosen [5], 
Ralston [8], and Fehlberg [12], 

Ralston [8] gives an optimum method (minimum TEB). 


0 

.4 

.45573725 

.4 

.29697761 

.15875964 


Ralston Optimum 
A. 21 

1 

.21810040 

-3.05096516 

3.83286476 



.17476028 

- .55148066 

1.20553560 

.17118478 


In King [11J, two optimum methods are derived, both whose TEB's 
are larger than Ralston's optimum, eq. A. 21. The first of these be- 

jy 

comes fifth order, and the second, sixth order when — is independent 
of Y. i 
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optimum, eq. A. 21. 


0 

2/5 

3/5 

2/5 

-3/20 

3/4 

F 

Kuntzmann 
A. 24 

1 

19/44 

-15/44 

40/44 



55/360 

125/360 

125/360 55/360 



Hull and Johnston [13] point out that, c_ « .35 anc i c ~ .45 

^ J 

(solve deriving equations for remaining coefficients) lead to a 
minimum TEB. As previously mentioned in A. 1-2, these bounds are 
usually larger than the actual error. 

Lapidus and Seinfeld [ 6 ] list four methods; one of these, 
the Gill form, is optimized for round-off errors (has a larger 
TEB than Ralston’s optimum). 




0 

1/3 

2/3 

1 



1/3 

-1/3 

1 


1 

-1 1 


1/8 3/8 3/8 1/8 


Method (p,V) 


Ralston A. 21 (4,4) 
King A. 22 (4,4) 
King A. 23 (4,4) 
Gill A. 25 (4,4) 
Classic A. 26 (4,4) 
Kutta A. 27 (4,4) 


Classic 
A. 26 


Kutta 
A. 27 


Error Bound (TEB) 
.0546 ML 4 h 5 
.0944 ML 4 h 5 
.1218 ML 4 h 5 
.0882 ML 4 h 5 
.1014 ML 4 h 5 
.0991 ML 4 h 5 


Tests made by the author showed that the Gill and Classic forms 
aie equal on accuracy while the Ralston form is superior to both these. 
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Hull [19] gives a good discussion on the problem of optimizing 
Runge-Kutta and predictor-corrector type methods. 

The ambivalence of the minimum TEB measure can be pointed out; 
according to Luther and Sierra [14], using their own minimum TEB 
measure, the Kutta form, eq. A. 27 is optimum for truncation error. 

Blum [15] claims to wipe out the optimum round-off error ad- 
vantage of the Gill form, eq. A. 25, by modifying the arithmetic 
(computer) sequence of the Kutta form, eq. A. 27, making it comparable 
to the Gill form. Fyfe [16] extends this procedure to any RKE(4,4) 
(though these modifications require extra programming effort). 

Lawson [17] derives a RKE(4,4) with an extended region of sta- 
bility; for integrating ordinary differential equations with large 
Lipschitz constants. This is comparable to the other RKE(4,4)'s 
on accuracy but can take larger step-sizes. Gates [26] gives a 
RKE(4,5) . 

Fehlberg [12] lists two methods with his usual feature of 
automatic step size control (see section A. 1-2). Also eqs. A. 16 
and A. 17, Fehlberg third order methods, have fourth order methods 


coupled. 

these can 

be used 

separately 

as fourth 

orders. 

0 





Fehlberg 

2/9 

1/9 




A. 28 

1/3 

1/12 

1/4 




3/4 

69/128 

-243/128 

135/64 



p=4 1 

-17/12 

27/4 

-27/5 

16/15 


P“5 5/6 

65/432 

-5/16 

13/16 

4/27 

5/144 


1/9 

0 

9/20 

16/45 

1/12 p=4 


47/450 

0 

12/25 

32/225 

1/30 6/25 p=5 
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worst. On time eq. A. 29 was the besf while eq. A. 30 was the worst. 
All these three methods, by the same test were faster and more accu- 
rate than the Kutta form, eq. A, 27. The time savings of these three 
methods is a result of their automatic step size control feature. 

Sarafyan [37] shows that the Classic form, eq. A. 26, ha's em- 
bedded within it first, second and third order methods. This is 
naturally expected of a fourth order method. This embedding prin- 
ciple of Sarafyan can be used to monitor step size. Using eq. A. 26: 

+ = Y n ^2 Y n +1 ava il a frl e by the complete 

calculation of the method. So here (second order accuracy 1/n) 

is available at no extra cost. Comparing Y^j and gives step 

[3] 

size control. Y n+1 is^not directly available, it is available at 

3 

•^h and not at h; so to use it for step size control requires extra 
effort. 

This embedding principle can be extended to any order method, 
and is a built in advantage of the Runge-Kutta methods. 

Various authors have studied the errors involved in fourth 
order Runge-Kutta methods and developed embedding methods to check 
these errors and hence institute step size control; regions of sta- 
bility have also been studied. These authors are: Merson [29], 

Scraton [20], England [21], Ceschino and Kuntzmann [32], Collatz [33], 
Lotkin [34], Chai [35], Sarafyan [36,37], Karim [38], Warten [29], 
Shampaine and Watts [40], Christiansen [41], and O’Regan [42]. 

Henrici [43] shows how to control round-off error. 

Computer programs using fourth order processes with automatic 
step size control are listed in Basnett [69] (uses the England [31] 
method), and in Jones [ 70 ] (uses the Classic, eq. A. 26, method). 
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The author, using the same principle outlined in section A. 1-3, 
has developed two new classes of RKE(4,4)'s: 

i- Back-step RKE(4,4) : 

0 Optimum 


-.1 

- .1 




A. 31 

.9 

4.65 

-3.75 




1 

7.44 

-6.24 

- .2 




1.06 

- .61 

1.11 

-.56 



Error Bound (TEB) < .22254 ML 4 . If for this method c 3 - .1 

and c 2 » .9 was used, the bound would be < .50641 ML 4 , more than 

twice that of eq. A. 31. This happens in every case, for the co- 

®£ficients of a back step method j so it is better to have c 2 as 

the negative back step coefficient rather than Cy One function 

evaluation per step is saved here as the k, of (T , Y ) can be used 

n n 

as the k 2 of ( T n+3 » ^ n+ ^)* This method is not self-starting, and 
needs a regular fourth order method with c 3 = .9 to start it. 


0 


-1/2 

-1/2 

1/2 

3/4 -1/4 

1 

- 2 12 


1/6 0 2/3 1/6 


Error Bound (TEB) < .26181 ML 4 . This method is close to the 
optimum, has simple coefficients, and can use either the Classic 
or the Gill form as a starting method. 
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Error Bound (IB, < . 23176 Hete ^ as ^ ^ ^ ^ 

Baok Step methods interchanging values of c 2 and c 3 results in 

a bound nearly twice as Urge, i.e., < .43719 ML 4 for = 2 . 1.01 and 

c 3 - For this method the k of 7 t v ^ , 

k 3 of (T n * Y n > can be used as the k 2 

° £ <In+1 ’ i ”+l > ' 0,16 m “ Jor advantage over the Back Step methods is 

that the Front Step methods are self starting. The optimum Front 

Step eq. A. 33 has a bound not much larger than the optimum Back Step 
sq. a. 31; so the Front Step methods should have great potential. 

As there are only two free parameters in the fourth order case, 
a back and front step (V - 2) fourth order method Is not possible; ’ 
but as the fifth order case has five free parameters a back and 

front step fifth order is possible. The author is working on such 
fifth and higher order methods. 

On comparing error bounds with the other fourth order methods 

ft is seen that the bounds of the Front and Back step methods are 

much larger; this is to be exDeet-AH - 

expected, and is the price paid for 

achieving V = 3 in a fourth order method. 


A. 1-5 RKE (5. 


The RKE(5,6)'s are claimed by many exponents 
field, to be the best compromise between accuracy 
ficiency (computer time). 


» in the Runge-Kutta 
and computational efr- 
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The deriving equations are presented in Butcher [3] , Fehlberg [4] , 
Luther and Konen [20], Luther [21], Konen and Luther [22], These also 
give general solutions to the deriving equations. The last three 
references explore the complete range of solutions possible, in terms 
of two and one parameter families. These lead to various Newton- 
Cotes, Gauss, Radan and Lobatto family of formulae. Cassity [23], 
Cassity [24] and Lawson [25] also state and solve the deriving equa- 
tions generally in terms of various free parameters. 

Rosen [4] gives deriving equations using the Quadrative approach. 

Lapidus and Seinfeld [6] also list a number of RKE (5,6) in- 
cluding the Kutta form corrected by Nystom. 


0 







1/3 

1/3 





Nystrom 

2/5 

4/25 

6/25 




A. 34 

1 

1/4 

-12/4 

f 

15/4 




2/3 

6/81 

90/81 

-50/81 

8/81 



4/5 

6/75 

36/75 

10/75 

8/75 

0 



23/192 

0 

125/192 

0 

-81/192 

125/192 


I . RKE (5 ,6) of Newton- Cotes Quadrative Family : . 


0 







1 

1 





Luther [21] 

1 

1/2 

1/2 




A. 35 

1/4 

14/64 

5/ 64 

-3/64 




1/2 

-12/96 

-12/96 

8/96 

64/96 



3/4 

0 

-9/64 

5/64 

16/64 

36/64 



7/90 

0 

7/90 

32/90 

12/90 

32/90 
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0 I 


1 

1 





Luther [21] 

1/2 

3/8 

1/8 




A. 36 

1 

-1/2 

-1/2 

4/2 




1/2 

4/64 

-5/64 

20/64 

-3/64 



3/4 

12/64 

9/64 

-12/64 

7/64 

32/64 



7/90 

0 

12/90 

7/90 

32/90 

32/90 



0 


1/4 

1/4 




Butcher [18] 

1/4 

1/8 

1/8 



A. 38 

1/2 

0 

-1/2 

1 



3/4 

3/16 

0 

0 

9/16 


1 

-3/7 

2/7 

12/7 

-12/7 

8/7 


7/90 

0 

32/90 

12/90 

32/90 7/90 


This Butcher form, eq. A. 38, Is recommended by Lapidus and 
Seinfeld [6] (after extensive tests) As the best RKE (p,v) to use. 
Also Sarafyan [44] uses this form for an ingenious error analysis, 




1-20 

hence this form is highly recommended* The author’s test show this 
method to be the best RKE(5,6) on accuracy* 



This Butcher form is of the back step type, discussed in section 
A* 1-3* This allows in all steps after the first, to use from 

( T n * Y n ) for the k 2 of ^ T n+l’ Y n+l^* He nce, overall, this method 
requires V = 5 instead of V = '6. Another such method is (not of the 
Newton-Cotes family) : 
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Here c. = -1/5 c- = 4/5 so allows the use of k_ from (T ,Y ) for 

5 n n 

the of ( T n +i» Y n+ ^)- T He larger number of zeros present in eq. 

A. 40 as compared to eq. A. 39 should make eq. A. 40 more computationally 
efficient than eq. A. 39. 
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The extra programming effort required to save the k's from a 
previous step, for use at a present step, is trivial. 

Another Butcher [18] method, which can be used to start eq. A. 40, 
is: 


0 


1/5 

1/5 




Butcher [18] 

2/5 

0 

2/5 



A. 41 

1/3 

7/36 

0 

5/36 



4/5 

0 

0 

4/5 

0 


1 

1/4 

0 

-35/4 

54/7 

25/14 


5/48 

0 

0 

27/56 

125/336 1/24 


Lawson [25] claims a form, similar to his RKE(4,4). This method 
has an extended region of stability: 
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The author’s tests confirmed that this Lawson form, eq. A. 42, 
does have a larger region of stability than other RKE(5,6)’s, and so 
can take larger step sizes. On accuracy this form is as good as the 
Bulcher form, eq. A. 38. These two are the best RKE(5,6)’s on 


accuracy. 
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_I1 . RKB(5,6) of the Gauss Quadrative family : 
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HI » RKE(5, 6) of the Radau Quadrative f a m ily : 
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The deriving equations for the Radau family of RKE(5,6)'s are 
given and solved in terms of free parameters in Luther [21] and 
Konen and Luther [22], The above equation is one such solution. 
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IV. RKE(5,6) of the Lobatto Quadrative family : 
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Konen and Luther [22] solve the deriving equations of the 
Lobatto family in terms of free parameters. The above two equations 
are two such solutions. 

% 

The advantage of the Newton-Cotes, Gauss, Radau and Lobatto 
-families is that when ^ is independent of Y, they become sixth 
order quadrative formulae of chese families. 
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Luther and Konen [20] point out that, on basis of computational 

! 

efficiency, due to the number of zeros present, the Gauss and Radau 
forms are better than the Lobatto forms, which in turn is better 
than the Newton Cotes form. While *a round-off error minimization; 
following the principle outlined in Gates [26], shows that on 
accuracy (based on round-off error) , the Lobatto form is better than 
the other three. 

Shanks [27] developed a class of RKE(5,5)’s rather than RKE(5,6)'s. 
This was done by solving the non-linear deriving equation approxi- 
mately, not exactly, resulting in one less function evaluation. The 
theory for reducing V for RKE's is also given here for p < 7, making 
it possible to derive more efficient formulae as compared to the con- 
ventional approach. However, it must be pointed out that (in the 
opinion of the author) methods derived this way are not exa.ctly p 
order, though they are greater than p-1 order. This is because the 
truncation error, in such methods has' a larger bound. The author's 
tests confirmed that Shank's RKE(5,5)'s are faster, though less 
accurate than, RKE(5,6)’s. 
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rs 

Shanks suggests that 1200 Nh > 1 be chosen to give a valid 
fifth order method. The author tested values of N from 100,000 to 
.0001 and found that N - 10 gave discrepencies in the coefficients, 
while N < 10 does give a valid method. Shanks suggests N = 9. 

The author’s test with various values of N showed that N = 5 gave 
the most accurate method. 

The author modified the Shanks form to give a method of four 
instead of five stages. This was done by using N = 100, resulting 

in c 2 = 1qq *q qq • As c 2 is much sraaller than c 3» c 4 or c 5 ; c 2 = 0 
was used. This makes k 2 = ^ and hence one function evaluation is 

saved. 



The author’s tests showed that this method is faster though 
less accurate than the other Shanks forms. 

Fehlberg [4] develops methods of order p, coupled with a p+1 
order method for automatic step-size control, based on minimizing 
the TEB for the p order method as discussed in section A. 1-2. 



Sarafyan [28] develops similar type methods with step size con- 
trol, based on the embedding principle outlined in section A. 1-4, 

But Sarafyan prefers to control his p order methods by embedded 

p-i, i = 1,2,3 , order methods. This requires no extra function 

evaluations, so is faster than Fehlberg* s methods. But as Fehlberg 
uses a p order method controlled by a p+1 order method, his error 
estimate will be better than Sarafyan; and hence a larger step size 
is possible with Fehlberg T s methods, at the cost of extra function 
evaluations. However the Sarafyan methods can obviously be used in 
the same way as the Fehlberg methods, by using the p order method to 
control a p-1 order solution. 



The Sarafyan and Fehlberg step size control principles are con- 
venient for use and are extendible to any order method, 

A Sarafyan fifth order with an embedded fourth order, is eq* A, 30 
listed in sections A. 1-4* Sarafyan [28] gives the deriving equations 

of such methods and solves them generally, giving six such embedded 
methods „ 
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15 
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All these embedding forms eq. A. 50 to eq. A. 55 are comparable 
to other conventional fifth order methods on accuracy but their 
step size control makes them faster* 

Sarafyan [44] employs an ingenious step size control with: 
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This form is gotten from the Butcher form, eq. A, 38, by multi- 
plying all coefficients by a factor of 2. Sarafyan shows that in 

going from (T , Y ) to (T - , Y ); using eq. A. 38 with h, Y 

n n n+1 n+1 n+1 

is available to a fourth order accuracy; at the same time, using 

eq. A. 50 with 2h, Y j available to a fifth order accuracy. 

Net result - overall V=3 is achieved although the RKE(5,6) of 

Butcher is used; and by comparing fourth and fifth order Y *s 

n+1 

step size can be controlled. 

In the opinion of the author, this Sarafyan method should be 
more accurate than eq. A. 50 to eq. A. 55, because the most accurate 
fifth order, eq. A, 38, is used here. 

A. 1-6 RKE (6 , 7 ) : 

The deriving equations are presented in Butcher [3], Fehlberg [4], 
and Rosen [5] . The first two also give the general solution of these 
equations, in terms of various free parameters. 
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In Butcher [3] are listed four RKE(6,7)'s 
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A back 

step and 

front step type method 

(see sections A. 1-3 and 


A. 1-4) 

i is also 

listed. 
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In eq. A. 58 k_ from (T , Y ) can be used for the k r of (T 

3 n’ n 5 n+1’ 

Y , 1 ); and the k, of (T , Y ) can be used as the k. of (T , . Y ,,). 
n-rl o n n 4 n+1 n+1 

So overall V=4 is achieved, hence this method is highly recommended. 




Luther [45] derives an ’optimum' Lobatto form: 
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All these quadrative forms eqs. A. 57 to eq. A. 61 become seventh 
dY 

order when ^ is independent of Y. 

Lawson [46 J lists a form with an extended order of stability 

(similar analysis to his RKE(4,4) and RKE(5,6), but admits that its 

accuracy improvement is marginal compared to other sixth order 

formulae, in fact worse for some non-linear problems, 

Fehlberg [4] lists a RKE [ (6 , 8) ; (7 , 10) ] , with the usual advantage 

of Fehlberg methods (see section A. 1-2), namely automatic step size 

control and minimum TEB. Sarafyan [47] lists four RK(6,8) with their 

embedding advantages (similar to his fifth order methods, see section 

A, 1-5). Huta's RKE(6,8) is included in this reference. Some, rather 

obscure, advantages are claimed for these so called improved sixth 

order methods. One definite advantage is that all four are of the 

dY 

Newton-Cotes family and when is independent of Y they become 
eighth order in accuracy. 
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V T n + h) ■ Y n + s£o [41(k 1 +k 8 ) + 216 <W + 27(k 4 +k 6 ) + 272k 5 ] 

where 

k x “ hf <W Huta 

k 2 - hf(T n+ ^, Y n+ ik x ) A . 63 

k 3 - hf( T * + |h. v n+ -L <k 2 + 3k 3 )) 

k 4 - hf (T n + 3k, Y n + | (k x - 3k 2 + 4k 3 )) 

k 5 - hf (T n + jh. Y n + A (-5k x + 27k 2 - 24k 3 + 6k 4 )) 

k 6 * hf(T n + I"* Y „ + ? (2Zlk i ’ 981k 2 + ®« k 3 " W2k 4 + k 5 )) 

k 7 " hf(T n + f* 1 * Y n + TS (_183k l + 678k 2 ' 472k 3 * 66k^ + 80k 5 3 kf> )) 

k g = hf(T n + h, Y n +I <716k 1 - 2079k 2 + 1002k 3 + 834k 4 - 454k,. - 9 kft + 72k ? )K 
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V T n + h) = Y n + m t 41 < k i +k 8 ) + 216 < k 3 +V + 27(k 4 + V + 272k 5 ] 

where 

k * hf (T Y ) 

± n n 

k 2 = + 9^» + k i^ Sarafyan 

k 3 = hf <T n + 6* 1 * Y n + < k l + 3k 2 » A.64 

k 4 - hf <T a + ^h, Y n + j (k x - 3k 2 + 4k 3 » 

k 5 = hf (T n + ¥• Y n + | (k l + 3 V> 

k 6 " hf (T n + ¥• Y n + I (-4k l " 21k 2 + 46k 3 " 29k 4 + 10k 5 >> 

k 7 = M (T n + Y n + Jl (_8k l + " k 2 " 84k 3 +44k 4 +9k 6 >) 

k 8 “ M (T n + h » Y a + J? < 107k 1 - 243k 2 + 354k 4 - 172k 5 -36k 6 + 72k ? )) 

k. = hf (T ,Y ) 

1 n n 

k 2 ” ^ ^ T n + "i* 1 * Y n + Sarafyan 

k 3 * hf(I n + £>■ Y „ + -k <k l + 3k 2» A.65 

k 4 ■ “ <T n Y n + i (k l - 3k 2 + 4k 3» 

k 5 ' ht(T n + ¥• \ * ¥ <k l + 3k «» 
k 6 ° hf(T n + | h - \ + - ?kj + -J (Wk 3 + k 4 » 

k 7 - hf <T n + |h, Y n + i (-SSkj + 99k 2 + 96k 3 - 180k 4 + 104k,. + 9k,.)) 

k 8 * h£<T n + h ’ Y n + ^ (287k l - 243k 2 ' 540k 3 + 894k 4 - 352k, - 36k 6 

+ Y2k ? )) 
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k l " hf <W 

k 2 = hf (T n + Y n + ? k l ) Sarafyan 

k 3 = hf (T n + |h, Y n + ± (k x + 3k £ )) A. 66 

k 4 - h£(T n + |h, Y n + \ (k L - 3k., + 4k 3 » 
k 5 - h £( T n+ ih, Y n + i (k 3 + 3k 4 » 
k = hf (T + -|h, Y + -i (17k. - 63k + 51k_ + k,)) 

o TiJn? 1 z j j 

k = hf (T + -fh, Y + -jy (-22k. + 33k„ + 30k- - 58k. + 34k. + 3k,)) 

7 n 6 n 24 1 2 3 4 5 6 

k Q = hf(T + h, Y (281k. - 243k„ - 522k_ + 876k. - 346k, 

o n n oz l z J h d 

-36k 6 + 72k 7 )) 


Shanks [27] developed an almost (6,6) method on the same 


lines as his RKE(5,5), eq* A*47* 
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A. 1-7 RKB(7,10) : 

Butcher [3] , Fehlberg [4] , and Rosen [5] list the deriving 
equations. 

Sarafyan [48] gives four such formulae, of his usual embedded 
types (see section A. 1-4), with built in step size control. His 
best form on accuracy is: 

Y T n +h > = V k l 

Y 2 ( T n +h) = Y q + ■|(-55k 1 +63k 2 ) 

Y 7 (T n +h) = V 17580 1751 (k l +k 10 )+357 7 <V k 9 )+1323 (k 5 +k 8 )+2989(k 6 +k 7 ) ] 

k = hf (T » Y ) 
i n n 

k 2 = k ^^ Y n + 1>3^* Y n + ^ k l^ Sarafyan 

k 3 = hf (T n + jjh, Y n + ^(k 1 +3k 2 )) A. 68 

k 4 " hf ' < T n + ¥> Y n + H (k l +3k 3» 

k 5 = hf(T n + |h, Y n + y(k r 3k 3 +4k 4 )) 

k 6 = hf ^ T n + T* 1 ’ Y n + 28^ k l +3k 5^ 

k ? = hf (T n + jh, Y q + -^^(-101^+651^-477^+449^)) 

kg = hf(T n +|h, Y n + 6577 ' 2 (-1881k 1 -783k 3 +10352k 4 -3414k 5 +5122k 7 )) 

k 9 = hf (T n + |h, Y n + 222—^(683663^+430650^-2032615^+2208930^ 

+385270k 6 -740735k ? +970137k g ) ) 

k 10 = hf ^ T n + h,Y n + 19601100 ^~ 12175/ * 21k l" 11236050k 3 +62891 ^ 30k 4 

~43488585k 5 -9947l40k 6 +5l099720k 7 -30879954kg+13337100k 9 )). 
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W h) - V k i 

y 2 (r o +h) - Y n + l(- k l«k 2 ) 

Y « ( V h) -vl (t i H W 

Y 7 (T tt +h) - Y n + “[41(k 1 +k 10 )+216(k 5 +k g )+27(k 6 +k 8 )+272k 7 ] 

k = hf(T ,Y ) 

1 n n 

k„ = hf(T + -rh, Y +%c-) Sarafyan 

z n j n J x 

k 3 = hf(T + -h, Y n + -|[k 1 +3k 2 ]) A. 69 

k 4 = hf(T n +h, Y n + |[k r 3k 2 +4k 3 ] ) 

k 5 ’ “ <T n + \ + 6*8 C83k 1 +32k 3 -7fc,. ] ) 

k 6 " hf <T n + I 11 ' Y n + ^[-3k ; -4k 3 +k 4 +24k 5 ]) 

k ? = hf(T n + -|h, Y n + 5^8 [-290^-524^+145^+1908^+1305^.] ) 

k. = hf(T + |h, Y + T 7^[292k.+108k,+13k,-318k c +753k,+106kJ) 

8 n 3 n 1431 L .1 3 4 5 6 7 J 

kg = hf (T n + |h, Y n + -£g|gg[ 14042^+11012^-4477^+5724^-6903^ 

+6360k_+31482k ] ) 
f 8 

i J 

k = hf (T + h, Y + T ^ r [-2Q49k.-1836k.+839k / +5724k -4692k, 

10 n n 4346 L 1 3 4 5 6 

+12084k 7 -9540k 8 +3816kg] ) . 
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W h > ■ V n^l 

VV h) ’ Y „ + |(-2Sk l +27k 2 ) 

[_VV h) * Y n + 84ol 4l < k x +k 1 o> +21 ‘( k s +k 9 > + 2?0‘ 6 +k 8 )+272k 7 i 
V T n + • Y „ ♦ & 

V T n + * Y „ + 1!<-V3 k 2> 

V T n + lh) . v n + ^(k l+ «k 3+ k 4 ) 


where 


Sarafyan 


A. 70 


and 


k i ’ hf <W 

k 2 - hf < T „ ♦ T? h * Y n ♦ ^ k !> 
k 3 • hf < T „ + Y n + ^(k 1+ Jk 2 )) 

k , - l»f( T „ + ih. Y n ♦ i(k r 3k 2 +4k 3 » 
k 5 " hf < T „ + l"’ \ + T4< k l +3k 4» 

k 6 “ hf <T n + 3*** Y n + ^(l«k 1 -420k 3 +411k 4 -128kj)) 

k 7 * l,f(T n + I* 1 * Y n + T^6{-2 4 9 4k 1 -H>876k J -5013k 4 +636k s +843k 6 >) 

k fl - hf<T n +|h. V n + T ^(-2456k 1 +7548k 3 -7977k 4+ 3074k 5 -189k ( . + I06k 7 )) 

k 9 " hf(T „ + f ) ’> Y n + 7^2(^946k I -118428k 3 +89043k 4 -7844k 5 -4431k 6 
-424k 7 +3498k fl )) 

k 10 * hf(T n * h - Y „ + 4^6<- 67 ”? k 1 +169308k 3 -106S33k 4 -3180k 5 +4272k 6 
+1378Ok 7 -954Ok 0 +3816k 9 )) , 


k 6 " hf(T n * 3* 1 * Y n + «(323k 1 -942k 3 +915k 4 -278k 5 )> 

k ? - hf (T n + ih. Y n ♦ X~ (-2070k x +5804k 3 -3741k 4 +212k 5 +843k 6 )) 

k 8 “ h£(T n + !*»•, Y „ + T558(- 244 90 k 1 +72132k 3 -70I25k 4 +23108k 5 
-201k 6 +848k 7 )> 

k 9 ’ hf<T n + &>' Y n + ■JJJJg(h0370k x -I51248k 3 +96438k 4 +8480kj-9747k 4 
+5247k g >) 

k 10 " hf(T n + h ’ Y n + 8^2 <-125084k l +297912k j- 1 50897k 4 .42824k 5 

. -*-1474Sk fe +26288k 7 -19O80k 8 +7632k 9 )) . 
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w h ' • V k i 


V 2 (T n+ t.) - V n -8k 1+ 9k 2 


W h) - Y n + l4Hk 1 +k 10 )+216(k J +k 9 )-»-27(k 6 +k 8 )+2?2k 7 ] 


Y i ( Vi h) ■ Vl k l 


VV i h > ’ Y n + i<- k l +3k 2> 
W l h) * Y „ + 4<V*W 


where 


k x * hf <VV 

k 2 ’ hf <V TP’ \ + !l k l> 

S * hf(T n + H h - Y a + ZS<V 3k 2» 

k 4 - hf<T « + J h * Y n +n <k l' 3k 2 +4k 3 )) 
k 5 - hf(T n + ±h, V n+ ^(k 1+ 4k 2+ k 4 )) 


Sarafyan 


A. 71 


and 

k 6 ■ hf CT n + 3 1 ** Y n + ^(6V 12k 3 +221 V 206k 5 )) 

k 7 “ hf(T n + 5 h * Y n + ilftMisk^akj-aaftk^is^)) 

k 8 “ h£(T n + I 1 *’ Y n + 4f7 (-l71k i +456k 3 - 83 98k 4 +8268k 5 -49k 6 +212k 7 )) 

k 9 - hf(T n + |h. Y n + Yj^g(17367k ;L -52248k 3 +36218k 4 +13X44k 5 -10188k 6 +5247k 8 )) 

k 10 ’ hf<T n + hf Y n + 2l73*~ 11407k i +35012k 3 -3163 l> t 4 -l 7 490k 5 +3979k 6 +6572k 7 -4770k a +190ak 9 )' 


k 6 * hf(T n + 3 *»* v n + |{2k 1 -4k 3 +109k 4 -104k 5 » 

k 7 - hf < T n + '| h » + 4j4<* 11 5k 1 +448k 3 -548k 4 +212k s +21Sk 6 )) 

k 8 “ h£(T n + !**• Y n + 477 ( ~ 171k l +456k 3 -12428k 4 + 12296k 3 -49k 6 +212k 7 >) 

k 9 “ h£ (T n + I* 1 * Y n + 3jL6 <5789k i- 17416k 3 +1 9 848k 4 -3392k s -3396k 6 +1749k 8 )) 

k 10 " h£(T n + h * Y n + 2D3 ( ' ll407k l +35012k 3- 21277 V 7844k 5 +3979k 6 +8 572k 7 -477Ok 8 +19O8k 9 )). 



Shanks [27] following the approach discussed in section A. 1-5 
lists a RKE(7,7) and a RKE(7,9). On a test system, as expected, 
the (7,9) was more accurate, while the (7,7) was faster. 
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A. 1-8 RKE(8,11) : 

Butcher [3] , Fehlberg [4] and Rosen [5] list the deriving equa- 
tions. 

Curtis [49] lists a RKE(8,11), listed here as eq. A. 73. In 
eq. A. 75 the 'b^'s 1 are the coefficients w^s used in this study. 
Curtis leaves C 2 > w^(bg), and w^O)^) free; and recommends - 
c 3 , Wg(b 9 ) = 1/5, and Wj^O^q) = 13/80. 
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Fehlberg [4] lists a RKE[ (8,15) ; (9,17)] with his usual advantage 

of automatic step-size control and minimum TEB (see section A. 1-2). 

Fourty digit arithmetic is used for the coefficients of this method. 

This is listed here as eq. A. 76, where the 's.'s*, B. .’s', and ' c.'s' 

i ’ ij i 

are the same as the coefficients, c^, a„ and w^ respectively, used 
in this study. Also in eq. A. 76 i and j take values from 0 to 16 while 
the formulation used in this study avoids i,j=0. Hence (3^ of eq. 

A. 76 is the usual a^^, is and so on. Naturally a^(c^) * 0 in, 
eq. A. 76. A formula for calculating the truncation error is also 


listed. 
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», = 
<*2 12 
Qfa 12 
cr 4 = 
Off. = 
Qt = 
Ob - 
Oa = 

Og = 

«10 = 
Qf U = 
0\3 = 

O' 13 = 

0-14 = 
Ofjs = 
<*16 " 

Bio = 
3a = 
3a — 
3a = 

333 ~ 

340 = 

343 “ 
343 = ‘ 

3 50 = 

3^ = 
3« = 
060 = 
3a3 = 

3e* = 
Bee =■ 
3to =• 
3t4 - 
Pts - 

3 76 ” 

3a - 


0.4430 8940 3704 9818 3103 0994 0428 1370 
0.GG55 3410 0647 4727 4CG4 3991 0042 2050 
0.9983 0115 8471 2091 199G 098G 0903 3083 
0.3105 0000 0000 0000 0000 0000 0000 0000 
0.5054 4100 9481 6906 8626 5161 2673 7384 
0.1714 2857 1428 5714 2857 1428 f.714 2807 
0.8280 7142 8071 4280 7142 8071 4 280 7143 
0.6054 3966 1210 1156 2034 9037 6920 0086 
0.2487 8317 9680 6200 20G9 7222 740G 0771 
0.1090 0000 0000 0000 0000 0000 000(1 0000 
0.8910 0000 0000 0000 0000 0000 0000 0000 
0.3995 0000 0000 0000 0000 0000 0000 0000 
0.6005 0000 0000 0000 0000 0000 0000 0000 
1 
0 
1 

0.4436 8940 3764 9818 3103 5994 0428 1370 

0.1663 8352 6411 8681 8GG6 0997 7660 5514 

0.4991 5057 9235 6040 0998 2993 2981 G541 

0.2495 7528 9617 8022 7999 1496 6490 8271 

0.7487 2586 8853 4068 3997 4489 9472 4812 

0.2066 1891 1634 0060 2426 5567 1039 3185 

0.1770 7880 3779 8634 7040 3609 9728 8319 

•0.6819 7715 4138 6949 4G69 377 0 7081 5048 • 10 -1 

0.1092 7823 1526 6640 8227 9038 9092 6157 

0.4021 5962 6423 6799 5421 9905 6369 0087 • 10"* 

0.3921 4118 1 690 7898 0444 3923 3017 4325 

0.9889 9281 4091 646C 5304 8447 6543 4355 • 10” 1 

0.3513 8370 2279 6396 6951 2044 8735 6703 * 10~ 3 

0.1247 6099 9831 6001 6G21 5206 2587 2489 

0.5574 5546 8349 8979 9643 7429 0146 6348 • 10* 1 

0.3680 6865 2S62 4220 3724 1531 0108 0691 

0.2227 3897 4694.7600 7645 0240 2094 4166 * 10 +1 

0.1374 2908 2567 0291 0729 5656 9124 5744 • 10 +1 

0.2049 7390 0271 1160 3002 1593 5409 2206 • 10 41 

0.4546 7962 6413 4715 0077 3519 5060 3349 * 10" 1 


Fehlberg 
A. 76 
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Pi4 6 =-0.6557 0189 4497 41 04 5138 0008 7998 5251 
0W7 =-0.3908 G144 8804 39SG 3435 0255 2024 1310 
0146 = 0.2670 4G46 7128 5002 2936 5844 2327 1209 

0149 =-0.1038 3022 9913 8249 0865 7698 5850 7427 - 

0i4io = 0.1667 2327 3242 5867 1664 7273 4616 8501 ■ 
014U = 0.4955 1925 8553 1597 70G7 7329 6707 1441 

01412 = 0.1139 4001 1323 9706 3228 5867 3814 1784 • 

01413 ~ 0.5133 6696 4246 5861 3688 1990 9719 1534 ■ 

0iec = 0.1046 4847 340G 1481 0391 8730 0240 6755 • 
01SB ="0.6716 388G 8449 9028 2237 7784 4617 8020 • 
0is> = 0.8182 8702 1894 2502 1265 3300 6524 8999 • 
0isio =-0.4264 0342 8644 8334 7277 1421 3608 7561 • 
015U = 0.2800 9029 4741 6893 6545 97G3 3)15 3703 • 
0isi2 =*0.8783 5333 87G2 3867 6639 0578 1314 5633 • 

0isi3 = 0.1025 4505 1108 2555 8084 2177 6966 4009 • 

0iso =-0.1353 6550 7861 7406 7080 4421 6888 9966 • 

0165 =-0.183 9 6103 1 44 8 4 82 7 037 5 044 1 9898 8231 

Pies =-0.6557 0189 4497 4164 5138 0066 7998 5251 
01*77 =“0.3908 6144 8S04 3986 3435 0255 2024 1310 

0166 = 0.2746 6285 5812 9992 5758 9622 0773 2989 

0i® =-0.1046 4851 7535 7191 5887 0351 8857 2676 • 

01610 = 0.1671 4967 6671 2315 5012 0044 8830 6588 • 

01611 = 0.4952 3916 8258 4180 8131 1869 9074 0287 

01612 = 0.1148 1836 4662 7330 1905 2257 9595 4930 • 

01613 = 0.4108 2191 3138 3305 5603 9813 2752 7525 • 
01615 = 1 


10 +l 

10 +l 

10 +1 

10- 1 

10" 2 

10 -2 

10 -2 

10" 2 

10 ~" 

10- 2 

10" 1 

10 +1 


10 +1 

10 +l 

10 +1 

10" 1 


Co = 0.3225 6083 5002 1624 9913 6129 0096 0247 • 10" 1 

Cs = 0.2598 3725 2837 1540 3018 8870 2317 1963 

Cg = 0.9284 7805 9965 7702 7788 0637 1430 2190 • 10~ 

C 10 = 0.1645 2339 5147 6434 2891 6477 3184 2800 

c u = 0.1766 5951 6378 6007 4367 0S42 9839 7547 

c^g = 0.2392 0102 3203 5275 9374 1089 3332 0941 

c l3 = 0.3948 4274 6042 0285 3746 7521 1882 9325 • 10 -2 

c l4 = 0.3072 6495 4758 6064 0406 3683 0552 2124 • 10 _l 


Fehlberg 
A* 76 contd. 
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Pes = 0.3254 2131 7015 8914 7114 6774 69G4 8853 
Bee 3 0.2847 66C0 1385 2790 8888 1824 2057 3687 

= 0.9783 7801 6759 7915 2435 8683 9727 1099 • 10 -2 
P*> = 0.6084 2071 0626 2205 7051 0941 4520 5182 • 10* 1 
Pas =-0.2118 45G5 7440 3700 752G 3252 7525 120G ♦ JO -1 
3* = 0.1959 655 7 266 1 7083 1957 464 4 9066 2983 

Psr^ =”0.4274 2640 3648 1760 3675 144 8 3534 2899 • 10 -2 

P»8 = 0.1743 4365 7363 1491 1955 3234 5255 8189 • 10 -1 

Pioo = 0.5405 9783 2989 3191 7355 7857 2411 1182 • lO* 1 

Pice ■ 0.1102 9325 5978 2392 6339 2831 2764 8228 

Pio? =“0.1256 5008 5200 7255 6414 1477 6378 2250 • 10 -2 

Pice - 0.3679 0043 4 775 8146 0136 3 84 0 4 33 6 6339 • 10~ 2 

Pi® =-0.5778 0542 7703 7207 3040 8406 2837 1866 • 10' 1 
Bno = 0.1273 2477 0686 6711 4546 6451 8179 9160 
Pu? = 0.1144 8895 0053 9310 5323 6588 7572 1817 
Pne * 0.2877 3020 7036 9799 2776 2022 0184 9198 
Piis = 0.5094 5379 4596 1136 3153 7358 8507 9465 
Pino =“0.1479 9682 2443 7257 5900 2421 4444 9640 
Piso =“0.3552 6793 8766 1674 0535 8485 4439 4333 • 10" 2 
Pis = 0.8162 989 6 0123 1891 9777 8194 2124 7030 • 10" 1 
Piste =-0.3860 7735 6356 9350 6490 5176 9434 3215 
Pia7 = 0.3086 2242 9246 0510 6450 4741 GG02 5206 • 10" 1 
Pi® =“0.5807 7254 5283 2060 2815 8293 7473 3518 * 10 -1 
Pi® = 0.3359 8659 3289 8497 1493 1434 5136 2322 
Piaio= 0.4106 6880 4019 4995 8613 5496 2278 6417 
Piai=“0. 1184 0245 9723 5598 5520 6331 5615 4536 • 10* 1 
8 iso =“0.1237 5357 9212 4514 3254 9790 9613 5669 * 10 +1 
Pios =“0.2443 0768 5513 5478 5358 7348 6136 6763 • 10* 2 
Pi® = 0.5477 9568 9327 7863 6050 4365 2899 1173 
P 137 =“0.4441 3863 5334 1324 6374 9398 9056 9346 ♦ 10 +1 
Pi® = 0.1001 3104 8137 1326 6094 7926 1785 1022 • 10* 2 
Br» =“0.1499 5773 1 020 5175 8447 1709 8507 3142 * 10+ 2 
Pi 3 io= 0.5894 6948 5232 1701 3620 8245 3965 1427 • 10 +1 
Biou= 0.1738 0377 5034 2898 4877 6168 5744 0542 * 10 +1 
Pi3i2= 0.2751 2330 6931 6673 0263 7580 2286 0276 • 10* 2 
Pi4o =“0.3526 0859 3883 3452 2700 5029 5887 5588 
0146 =“0.1839 6103 1448 4827 0375 0441 9898 8231 


Fehlberg 
A. 76 contd. 
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Increasing use of higher precision arithmetic is required, to 
realize the advantages of this order over the seventh order (slim 
advantage in accuracy) — paying the price of increased computational 
time. 

A, 1-9 RKE T s of Ninth and Higher Order 

The deriving equations for ninth and higher orders are listed, 
but not solved in Rosen [5] . 

Shanks is supposed to have developed (not yet published) ninth 
and tenth order RKE's. 

There is no real reason for pushing on for orders beyond the 
eighth (except fot special problems) . According to Fehlberg [4] 
accuracy hardly improves after the seventh order methods. 

Saiafyan is working on a computer formulation and solution, 
of the deriving equations leading to high order methods (see 
Sarafyan and Brown [7]); however the results are yet to be published. 

A. 2 Semi- Implicit Runge-Kutta Methods - RKSI (p,V) : 

It is easier to use the autonomous form of an ordinary differential 
equation system when dealing with these types of methods. 

dY 

^ = F(T,Y) A. 79 

The explicit dependence of T in F(T,Y) can be removed by making 

ft 4* * 

T the (n+l)th component of the n component vector Y_. This makes the 
system n+1 dimensional rather than n dimensional. 


let 
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dY 


n+1 


dT 


= 1 


dY 

' • dT " - - ) 


A. 80 


where: 


/»! \ 


dY 

dT 


and F(Y) = 



f a • 
r i 


•Y •••Y Y ,_) 
i n n+1' 


W'VVW 


VY 


•Y •••Y Y ) 
i n n+1 




The initial conditions would also be accordingly altered. 


i,e ” WV 


Tq. There is no need however to solve the n+1 


component equation of eq. A, 80 as its solution is always Y^ + ^ = T. 

When a Runge-Kutta method is used on this autonomous form; 

the increments coefficients of the abcissa, i.e., the c's, become 

a’s, i.e., the increment coefficients of the vector component Y^ + ^. 

In Chapter 2, section 2.2 the method of solving for each k. 

of the semi- implicit form was discussed; this involved a Jacobian 

inversion. Restating eq. 2.33 in matrix form to apply to a system 

of equations, with the a_'s involved in the Jacobian written as 

a ! . ' a . 
ij 

i-1 1 i-1 

k. « [I - h a. . F (Y + l a’ k.J • hF(Y + T a. . kj 
-i ii — y — n ij -J J ~ -n ij -y 


A. 81 


i = i )2 , V 


The RKSI type method used to solve eq. A. 80 is: 
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V 

Y = Y + I w.k. 
— n+1 —a . l— i 
i=l 


k, - W(Y, +■ l ay 1^) 


j-1 


A. 82 


To solve for each k^, the formulation of eq. A. 81 is used. 
The coefficient matrix is of the form: 




Allen [50] » Rosenbrock [51] and Haines [52] list the deriving 
equations. Allen [50] gives a thorough discussion of the deriving 
equations and the stability advantages of RKSI’s over RKE's; and also 
gives coefficients of the truncation errors. In general RKSI’s re- 
quire less stages, V per step as compared to RKE's, to achieve a 
p order method. The price paid by RKSI's, for achieving extended 
stability as compared to RKE; s , is the additional computational ef- 
fort required to compute Jacobians and their inverses; but this price 
is not ds steep as the one paid by RKI's, where iteration for k^'s 
is required. 

Lapadis and Seinfeld [6] recommend the use of RKSI only for 
those cases where stability is critical, for example, with still 
systems; because in general RKSI's are not significantly more accurate 
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than comparable order RKE's; and RKSl's usually require more compu- 
tational time as compared to RKE's. 

The following sections are titled EKSI(p) instead of RKSI(p,V). 


A. 2-1 RKSI(l) : 

Allen [50] states that no such case of interest exists. 


A. 2-2 RKSI (2) ; 

Allen [50] and Rosenbrock [51] list the deriving equations and 
solve them. Allen [50] gives a V=2 method while Rosenbrock gives a 
V=2 process. 


1/2 

1 


Allen 

0 

A. 83 



- 1 ) 
2 


0 



1 


0 


Rosenbrock 
A. 84 


Allen [50] however admits that his method, eq. A. 83 does not 
possess significant stability advantages over a second order RKE. 


A. 2-3 RKSI (3): 


The deriving equations are listed and solved in Allen [50] , 
Rosenbrock [51] and Haines [52] . 

Allen [50] lists three methods with V=2: 


0 

0 1/3 

-1/2 3/2 


1/3 


Allen 
A. 85 
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Calahan [53] also lists a V=2 method, which is the same as Allen's 
eq. A. 87. 

Haines [52] derives a V=4 method, the extra stages per step 
being used to give a truncation error estimate and hence allow step 
size control. This method consists of two equations, each with V=4 
stages. These two equations are compared to give a truncation error 
estimate. Both these equations have mostly common coefficients. This 
method is listed in two parts because each part can be used inde- 
pendently to give a third order method, or in combination to give 


an error estimate. 
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: ~1 

1 

1/2 

1 

1/2 

1 


1 

1/2 

1/2 

Haines 
A. 89 

2/99 

95/99 

2/99 

2/3 

0 

0 0 


19/9 

-43/18 

28/9 

-11/6 





1 

1 1 

1/2 1/2 1 

21/44 19/44 1/11 1/2 

10/3 -1 7/3 -11/3 

Error^ - (29/9) [ (Y with eq. A. 89) - (Y with eq. A. 90)] 

Haines [52] also compares his method, which combines eqs. A. 89 
and A. 90 to give step size control, with a RKE(4,4) which uses a simi- 
lar step size control. His results point out that his method is 
superior on basis of computing time and stability. 

A. 2-4 RKSI (4) : 

Allen [50] lists thirth fourth order methods with V=3. Three 
such methods are presented here, these three were judged in Allen [50] , 
and found to be superior to the rest of his methods on some points. 

4444779 Allen (minimum 

computational 
4444779 0 effort) 

A. 91 
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.28915127 



.57118895 

.28915127 


.53873135 

-.024693209 

.33215386 

.72862466 

4.21831023 

-3.94693489 


0 Allen (largest 

region of stability) 

0 0 

A. 92 


.86234511 

-1.28325229 1.23119918 0 

- .80332019 - .091154713 1.23119818 0 0 

- .32266573 - .85621119 2.17887692 


Allen (best for 
numerical stability) 

A. 93 


Butcher [2] has listed a RKSI(4,3); his method does not apply 
to the autonomous form of the differential equations, but to the usual 
form used throughout this study. (The method could easily be con- 
verted to apply to an autonomous form by simply removing the c^’s.) 


0 

0 



1/2 

1/4 

1/4 


1 

0 

1 

0 


1/6 

2/3 

1/6 


Butcher 
A. 94 


A. 3 Implicit Runge-Kutta Methods - RKI (p,V) : 
Restating the equations for a RKI(p,V): 


Y . = Y + y W.k. 

n+1 n . L . x i 


i=l 


A. 95 


V 


k. = hf(T + c.h, Y + 7 a ,k.) 

i n i n ij y 

The coefficient matrix for these types of methods is full; and 
hence each k^ has to be solved for iteratively at each integration 


step. 



1-56 


The use of these RKI’s is recommended for two situatipns: one 

where stability is critical and a self-starting process is desired; 
two when the differential equation system is such that the cost of 
computing f(T,Y) for a given value of T is high compared to the cost 
of repeating this computation with Y changed but T unchanged. In 
the second situation the iterative nature of an RKI would not be an 
objection; especially if the differential equation system is linear, 
then the iterations can be replaced by standard linear algebra 
techniques. 

Advantages: Stability characteristics better than comparable 

order RJCSI’s or RKE's. Always stable for any h. 
Higher p for a particular V as compared to RKSI’s 
and RKJE's. The deriving equations are easier to 
formulate than those of RKSI's and RKE’s. 

Disadvantages: For even moderately complex differential equa- 

tion systems the iterative solution of k.'s re- 
quires more computational effort than comparable 
order semi- implicit or explicit methods. 

The coefficient matrix will be presented in the following form: 



Most currently available implicit processes were derived by 
Butcher. The theory behind implicit methods, attainable order, etc.. 
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is dealt in Butcher [1] and Verner [54]. Deriving equations and 
implicit methods are listed in Butcher [2], Butcher [3] and Butcher 
[55]* Lapidus and Seinfeld [6] also list some of these Butcher 
methods * 

Implicit Runge-Kutta methods are based on quadrative formulae 
(The Taylor Series analysis can also be used to derive implicit 
methods but the quadrative approach is simpler) . Basically Implicit 
Runge-Kutta 1 s can be divided into three classes of methods: Gauss- 

Legendre quadrative forms, Radau quadrative forms, and Lobatto 
quadrative forms* 


Quadrative Form 

Abcissa points (c.'s) 
specified 1 

V 

p 

Gauss-Legendre 

All c. f s found as roots 
of a Legendre poly- 
nomial, no specified 

V 

2V 

Radau^ 

c^ = 0 specified 

V 

(2V-1) 

Radau^ 

= 1 specified 

V 

(2V-1) 

Lobatto 

= 0 and c^. = 1 
specified 

V 

(2V-2) 


when c^ = 0 

a ll a 12 ~ a 13 

when Cy = 1 

a lV = a 2V = a 3V 


IV 


0 


vv 


0 


Hence the Lobatto forms have the first row and last column of 
the coefficient matrix as zeros. Thus as = 0 and c^ - 1 were 
specified two k_/s, k^ and are available explicitly and so can 
be solved for by iteration* The Radau^. and Radau.^ forms have only 
one k^ available explicitly as only one c^ is specified. While the 
Gauss-Legendre forms have no k^ available explicitly as no c^ is 
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specified. Because of the number of abcissa points specified, the 
Gauss-Leg end re forms have the highest p for a given V, while Radau 
forms are next and Lobatto forms last. 

Butcher [2] also suggests a convergent iterative procedure 
applicable for calculating the implicit k/s of any RKI. 

If n is the iteration count for a particular k^ the convergent 
iterative procedure is: 


kf n) = hj(T + c h, Y + l a kf n) + £ a..kf n-1) ) 

1 n 1 n ijJ 3 

Butcher proves that this procedure is convergent. 


A. 3-1 RKI(p,V) l s of the Gauss-Legendre Quadrative forms : 

These kinds of RKI’s are listed in Butcher [2] from p = 2 to 
p * 10 (p = 2V) , an error analysis is also included. 


RKI (2,1) 
A. 96 


1 /3 . 

1 /3 


4 " 6 

1./3 1./3 

1 

2 + 6 4 + 6 

4 

1 

1 

2 

2 


RKI (4, 2) 
A. 97 



RU(6,3) 
A. 98 


%_ 2 

0 0 



, Vu 8 

r 30 9 


_ v^is ». yift 

lfi ' 3ft ” 30 

a JL 

0 38 ~ 24 

. yi» * 

+ 15 3ft 


S 4 

9 9 


ft 

IS 


RH(8,4) 

ui ^ ^ ^ A* 99 

m' + * + * + «*'-- 

«n* + m + m <*+*+ — "* *i 



33 




RKK10.5) 
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A. 3-2 RKI(p.V) 's of the Radau Quadrative forms : 

These types of RKI's are listed in Butcher [55], from p*l to 
p=5 (p = 2V - 1). Butcher [55] also describes an ingenious error 
analysis applicable to the Radau and Lobatto (see section A. 3-3) 
forms of RKI's. 

RKI(1,1) 

A. 101 



000 RKI (3,2) 

211 A. 102 

3 3 3 

1 1 

4 4 


0 

0 

0 

0 


o — v 7 o 

» + \Ai 

24 + \A\ 

108 - 73 v/0 

RKI(5,3) 

10 

74 

" )2i) 

«X) 

A. 103 

+ V f> 

o - v'V* 

108 + 7;i vTi 

24 — y/ti 

10 

7Ti 


120 



1 

10 + N'Vi 

H> — 



" :«• :ui 



RKI (3,2) 

A. 104 

RKI(5,3) 
A. 105 
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S -W -+-.W-W 






2*' 


1 

•*-— 0 
00 

0 


L _ it~y ? . _ !i+-y? __ . /TTwi , m A /i- 2^7 

[_ 120 ’ ^ 120 ■ ^ y 21 ’ \ 21 ’ 

1-hs/i , 7 + 5v A 7 , 2 + 5V7 , 2-6y/j , 

“■*—«« ? m —w * wr" 4 **’ 

N+ \/7 , 14 - Vi 22+17^7 , 22 — 17 \/l . 

** _ ' 37«o ^ _ :h«7" Vl!l - «• 5S~“- 

22 — 5\/7 , 22 + S>/7 , I - vT , I + y/l ,”| 

**-- 55— <*. w 


0 

0 

0 

a 

0 

0 

0 

0 


L-<* 

2 

Mi 

1 

** 84 

<*i — tiH-j-eV 

m*— o* 

M| 

1 

**“84 ** 

0 

RKI(12,7) 

2 

1 

2 

«/ 

32 

M*'-<V+M. 

M* + M,* 

, 1 

W 

U_ 

100 

4 t 

m'—W—m 

0 

0 

A. Ill 

HV 

2 



W '-S +W 

m»'+*V 

, 1 
*“5 


0 


i±<+ 

2 

«• 


u*+M*4W 

MT+<4% 

Ml+tfi-W 

1 

0 


1 

0 

«1» 

Ml* 

8^ 

25 

/ 

M|| 

Mil 

0 



1 

42 

U 

2m/ 

m 

625 

2m/ 

2m, 

42 



L.- 124 ~ 7v/T ^ + _ ,/l5+ 2VT6 /lB - 2 vT5 

L 1400 ’ 1400 ’ y 33 • y 33 ' 

_413 , , 413 W + WTS . M — 4\/l6 00 + 27V1S 

“*“ 303 “'' sift- • - ■— sm~' mi 

, W» - 27vTa , OW + lOOv/iS . <158 - lOOv/15 . 327(1 - IOOx/Is 

* wT~ m ’ m> "*• u55 — -• ^5 — - 

, 3270 + lOOv'lS ai-4y'i6 . 30 + 4vTs , 8 - vlfl 

** “ 27285 ' iiS ** W~ m - 

, 8 + vTis , Hoi-Wy'iS . an+73 V / i6 m-nv'Is 

-t - — 7„. ">■ is,-.. “ ‘ > “*• " — «..» — > ">* ^ "»• 


, iii + i7%/i5 , 


, 51+2^/16 , 61—2 v^TSl 

' "" 300 ’ 3U0 J ‘ 
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A. 4 Multistep Runge-Kutta Methods - RKMS(p.V) : 

One way of reducing the stages V of a p order Runge-Kutta is to 
use an implicit Runge-Kutta — the price paid would be the iterations 
required. Another way of reducing V for a given p order Runge-Kutta 
is to utilize the solutions of one or more previously evaluated steps. 
This result is a Multistep Runge-Kutta. The price paid here is that 
the inherent self -starting, easily institutable step size change, 
accuracy, and stability advantages of the single step Runge-Kutta’ s are 
lost. Another disadvantage (trivial) of a multistep method, compared 
to a single step method, is that extra storage is required to save 
the solutions of previously evaluated steps. The typical complex 
error analysis of single step Runge-Kutta’ s occurs in RKMS’s too. 

On the whole standard multistep predictor-corrector methods 
perform better than RKMS's; so there does not seem much point in 
using a RKMS instead of a RKE, if a multistep method is desired, 
a standard predictor 7 corrector is recommended. 

Based on the usual formulation of a Runge-Kutta method, applicable 
to a first order ODE, a RKMS(p,V) is defined by: 


S n V 


n+1 


“ Y n + l I "i 
“ 1 1 


s=l 1=1 


n-s+1 , n-s+1 
k i 


A. 112 


M. 


^n-s+l m . • /•£ 
i n-s+1 1 w i 


, n-s+1 , v . v 

+ h » Vs+l + l a ij 

j“i j 


n-s+1 , n-s+1. 
R ij } 


S — n+1 
n 


M ± £ V 


1 * 2 , 3 , * * 
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The superscripts in eq. A. 112 are used to identify w^, k^, 
and a^ at a particular time step. For example, 2^=1 makes eq. A. 112 
a single step Runge-Kutta; S^«=2 makes eq. A. 112 a two step Runge-Kutta, 
and the k^s of the (n-l)th time step are used together with the k^s 
of the nth time step, to arrive at the solution Y n+ ^. If S n =n+1; then 
the wolutions of all previously evaluated steps, including the initial 
conditions; are used to calculate Y ; Henceforth RKMS(p,V,S n ) will 
be used to designate these methods. The value of identifies the 
number of steps used in the RKMS. Bryne and Lambert 59 show that 
there is not much point in going beyond a two step (S^=2) Runge-Kutta, 
because both accuracy and stability of these RKMS's deteriorate as 
compared to single step Runge-Kutta' s. 

The rest of this section will discuss S =2, or two step RKMS's, 

n 

eq. A. 112 is reformulated in a more convenient form for s n *2. 


Y =« Y + J u.g. + 7 w.k. 

n+1 “ W 11 i-1 11 

*i ■ + V* Vl + £ Vi> 


A. 113 


k. = hf (T 


Mi 

+ Y , + “ijV 


M i SV 

Usually M^ = i-1 is used to give an explicit formulation. 

Rosen [62] gives the deriving equations for RKtyS(4,3,2) and RKMS 
(5,5,2) using the convenient quadrative approach (convenient as com- 
pared to the Taylor Series approach). Butcher [56,57] also gives 
the deriving equation, but the methods of these papers are really 
hybrid methods, not RKMS's, and so come under section A. 5. 
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Bryne and Lambert [59] present the deriving equations for two 
RKM's: a RKMS (3,2,2) and a RKMS(4,3,2). The general solution, error 

analysis and proof of convergence is also presented. Bryne [58] gives 
the parameters which make the above two methods optimum (minimum TEB 
based on the Ralston [8] criteria, see sections A. 1-3 and A. 1-4). 
Lapidus and Seinfeld [6] also list these two methods. 


RKMS (3,2,2) : 


n-1 

n 

o 

El 

H 

o 

It 

O 

b 2 -c 2 

C 2~ C 2 

d 21 =c 2 

a r c 2 


(18c 2 -5) 

u r 1_w i 

yf = — — , 

1 12c 2 

U 2 =_W 2 

W 2~ 12c 2 

c 2 free, c 2 ■ y gives 

an optimum method. 


A. 114 


RKMS (4,3,2): 


n-1 


o 

II 

c^=0 

b 2~ c 2 

C 2 =C 

b 3 =c 3 

C 3 =C 

d 21 =C 2 

a 21* 

d 31 =a 31 

a 3l“ 

d 32 =a 32 

a 32 = 

u 1 =l-w 1 

W i = ' 

u 2 =_w 2 

w 2 = ■ 

U 3~ _W 1 

w 3 = - 


n 


c 3 -a 


32 

[2c 3 (c 3 ~ c 2 ) ] 


A. 115 


i. i. 

[4+18c 2 c 3 -5(c 2 +c 3 )] 


12c2C 3 

[4-5c 3 ] 


2 ~ [12c 2 (c 2 -c 3 )] 

^ 5c 2~ 4] 
[12c 3 (c 2 ~c 3 )] 


J it J 

C 2 and are free, c 2 = *-*^ and 722779927 give an optimum method. 
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Lapidus and Seinfeld 6 also list the above two methods. 

These two methods, eq. A. 114 and A. 115, are not as accurate as 
comparable order RKE’s (though they are faster). 

Chat 60 develops an optimum RKMS(4,3,2) using an ingenious 
derivation. Chai notes that if an RKE(4,4) has W2=0 then the can 
be reformulated in a multistep form without changing the c^'s, a^'s, 
and w^’s of the RKE(4,4). The method is hence defined by the usual 
RKE(4,4) formulation but with 1^ given differently. 

k, = hf (T ,Y ) + h 2 c,f r (T.Y ) 

2 n n 2 n n 

A. 116 


+ 

2 ! 


<w * i; 


e 2 f 


""(T ,Y ) + 
n n 


where 


eq. 


e^ and e 2 are constants. 

By using a finite difference approximation for f'(T n , 
A. 116 Chai arrives at: 


Y ) in 
n 


k 2 = q^ + q^! + qjgj 


A. 117 


where: 


VY . 

>* 

1 

D 

n-1 

n 

n- 

. ’ . VY , 

V 

- 1 

w,g 

n-1 

i^l 

i . 


A. 118 


the g^'s being the k^'s of the previous step. 

Thus V**3 is achieved because no function evaluation is necessary 


for k 2 > 
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By comparing a Taylor series expansion of eq. A. 117, up to second 
order terms, with the first two terms of eq. A. 116, gives a set of 
conditions for the q's. 


q l + 


q 2 + q 3 


2 



A. 119 


Thus one q is arbitrary. 

An optimum RKMS (4,3,2) is hence formulated by Chal (minimum TEB) 
this uses the same c^s, a^'s an d w^s as the RKE(4,4) from which it 
is derived, hence no separate starting method is necessary. 


RKMS (4,3,2) : 
C l“° 

C 2 = C 2 q l 
c 3 = 1/2 

C 4“ 1 
- c 

a 

- 

1 


i ^ 26 

1 + T C 2 


42 

q 2 " " 5 C 2 


16 

q 3 " 5 2 


21 "2 

1 1 

2 ” 8c 


A. 120 


41 2c. 

1 

w-j^ = 1/6 


2 
- 1 


32 8c, 


a 42 ' “ 2c, 


w„ 


“43 z 
2 

W 3 " 3 


w. 


1/6 


c^ “ 1/2 or 1/4 is recommended. 


Chai’s tests show that this RKMS is superior on accuracy and time 
to the RKE with which it shares its coefficients. An error expression 
for corre-ting the solution, to yield a fifth order accuracy is also 
given. DUe to its simplicity of application and other advantages, 
this method is highly recommended. 

Gruttke [61] gives the deriving equations for a RKMS (5, 4, 2). 


These equations are solved to yield various methods which are 
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extensively 

tested. An optimum based 

determined. 


RKMS(5,4,2): 


n-1 

_n 

b 1 » 0 

C 1 = 0 

b 2 = c 2 

c 2 = C 2 

b 3 = C 3 

C 3 = C 3 

b 4 = c 4 

62 

c 4 = 85 

d 21 = c 2 

a 21 = C 2 

d 31 = a 31 

a 3l = c 3 " a 32 

d 32 = a 32 

31 - 40c 4 

a 32 ° 240 w 3 c 2 (c 3 - c 4 ) 

d 41 = a 41 

a 41 = C 4 - a 42 “ a 43 

d 42 = a 42 

1 

a 42 ~ [120 w 4 c 2 (c 4 - c 3 l 


d 43 ~ a 43 


U 1 = 1 " W J 


U 2 -W 2 


u 3 = " w 3 


- w 4 


[ 


A. 121 


<c 3 - c 4 )(31 - 60c ? > 
3(c 3 - c 2 ) 


(31 - 40c 4 ) 


■] 


31 - 60c 2 

a 43 = 360 w 4 c 3 (c 3 - c 2 ) 

w 1 = -| - ( w 2 + w 3 + w 4^ 

31 + 50c 3 c 4 - 40(c 3 + c 4 ) 
W 2 = 120 c 2 (c 2 - c 3 )(c 2 - c 4 ) 

31 + 50c 2 c 4 - 40 (c 2 + c 4 ) 
W 3 + 120 c 3 (c 3 - c 2 )(c 3 - c 4 ) 

31 + 50c 2 c 3 - 40 (c 2 + c 3 ) 

w, = 


4 120c 4 (c 4 - c 2 ) (c 4 - c 3 ) 

c 2 and c 3 free, c 2 = .23 and c 3 = .869 give an optimum method. 

Again, eq. A. 117 would be faster than comparable fifth order RKE’s 
but not quite as accurate. 
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There would be very little advantage in deriving an implicit or 
semi-implicit RKMS form (M i = V or M ± = i in eqs. A. 112 or A. 113), 

t 

because of the additional complexities involved. 

A. 5 Special Types or Hybrid Runge-Kutta Methods 

Two most widely used methods of numerically solving ODE's are 
the Runge-Kutta methods and the predictor-corrector methods. 

For complex systems of ODE's numerical solutions, using the usual 
Runge-Kutta's discussed so far, have one big disadvantage — function 
evaluations take up the major portion of the computational time in- 
volved. Multi-step predictor-corrector methods require considerably 
less function evaluations per step. So obviously reducing the number 
of function evaluations (V) of a Runge-Kutta, while retaining some 
or all the good points of a Runge-Kutta, would yield good methods. 

Some compromise between Runge-Kutta's and predictor-correctors is 
indicated. 

Single step Runge-Kutta's have certain definite advantages and 
disadvantages, and so do the predictor-corrector multi-step method. 


Characteristics of Comparable Order Methods 


Single Step Runge-Kutta 

Advantages : 

1) Self-starting 

2) Comparatively good stability 
characteristics because single 
step (so large h possible) 

3) Step size and order easily change- 
able at each step. 


Multi-step predictor-correctors 

Disadvantages : 

1) Not self-starting 

2) Comparatively bad stability 
characteristics because multi- 
step (h comparatively small) 

3) Step size and order change not 
easily instituted at each step. 
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Disadvantages: 

1) May require greater computa- 
tional time as comparatively 
large number of function evalua- 
tions required 

2) No easily computable and ac- 
curate error estimate available 

3) High order methods involve ex- 
tremely complicated deriving 
equations. 

Hybrid Runge-Kutta *s have been developed to combine the advantages 
of the Runge-Kutta and predictor-corrector methods, while trying to 
avoid their disadvantages. 

Rosser 63 developed a class of ingenious hybrid Runge-Kutta’ s 
with great potential. These methods favorably combine the advantageous 
features of Runge-Kutta* s and predictor-correctors. They may be de- 
scribed as implicit, multi-step, predictor-corrector Runge-Kutta' s. 

These methods are "Block" methods and proceed by blocks of N 
steps. Each block is completely independent of other blocks. The 
solutions from previously evaluated steps are only required from, 
and used within, a particular block. A large degree of flexibility 
is built into these block methods giving them three concrete advantages: 

1) Step size is constant within a block but can easily be changed from 
block to block 

2) Order of accuracy can be varied from block to block by changing 


Advantages : 

1) Computationally fast as 
comparatively small number 
of function evaluations re- 
quired (provided not too 
many iterations done) . 

2) Fairly easily computable 
and accurate error estimate 
usually available. 

3) High order methods easily 
deroved/ 
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N, the number of steps constituting a particular block; and by 
changing the number of iterations used in each block. 

3) Accurate estimates of errors incurred within a block are available 
and so give good pointers for the best values of h and N for the 
next block. 

These "Block” methods result in a considerable savings of function 
evaluations and are hence highly recommended. The Rosser "Block" 
Runge-Kutta' s can be defined using the following formulation: 

For a block of N steps: 


X 1 ~ *0 = X 2 ’ X 1 


x - x . = h > 0 
n n-1 


The subscript 0 denotes the starting or initial step, hence 
Xq, y^ are known (from initial values or the results of the 
previous block). 


y 0 -y 0 

V ■ 




Starting values 

for a block A. 122 


y m “ y 0 + mhf(x 0 ,y 0 ) 

y' r - f(x , y*) 

m mm 

r £ 0 m = 1,2,3, N 

The superscript r is the iterate count. 

' y 0 + h[ "l y 0 + a ll y i r + ' + a lM y m r + + W.* 1 


A. 123 
A. 124 


r+l , . r ' , r . 

y “ y n + hlw y. + a ,y, + 

•’m J Q L nrO ml- 7 ! 


+ a 


mm m 


r . . r n 

y + * • • + a y J 


mn'n 


A, 125 


y r+ ^ = y ft + h[w yl + a ,y, r + **• + a y r + • • • + a y r ] 
J n 7 0 L n-'O nl 7 l nnrm mrn 


r - 0 m * 1,2,3, * • *N 
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w, , w„ * * * w and a- .. , a, * • • a are constants defining a 
12 n 11 lm nn 

particular block method. 

The working of the block method is easily understood. Consider 
the r=0 case, i.e., the first iteration: y^ is calculated from 

t ^ O t 

eq. A. 125; this requires N function evaluations — y^°, y 2 > *’’’y n ° 

which are calculated from eqs. A. 122 and A. 124 (y^ is assumed known 

from a previous block and remains the same for any iteration; though 

yQ must be calculated for the first block from eq. A. 122 and hence 

the first block needs N+l function evaluations). Similarly y 2 » 

y^ . . • are calculated from eq. A. 125, these use the same y!°» 

y^ 0 ,,*** y^° used by y^. Thus a complete sweep for a block (one 

complete iteration of eq. A. 125) requires N function evaluations. For 

2 

r=l, i.e., the second iteration: y^ is calculated from eq. A, 125 

1 1 1 

in the same way as for the r=0 case: the values y^, y^> ■** y n are 

substituted in eq. A.,124 to yield yj 1 2 , y' 1 , *** y R , which are then 

2 2 2 

used in eq, A. 125. jy Jy • • • Yn are similarly calculated. 

Up to a point increasing r, the number of iterations, increases 
the order of accuracy of the solutions y^, * * • y n * For example, 

for a N=4 method (Milne method) with 0 < r < 4 each iteration gives 
y^ to an order of accuracy of 4+2, For r > 4 no improvement in or- 
der of accuracy is gotten. If r=3 is used a fifth order method re- 
sults and uses a total of 16 (4N) function evaluation for a block of 
four steps (N=4). Thus a fifth order Runge-Kutta is gotten which 
uses four function evaluations per step. If r=4 is used a sixth 
order Runge-Kutta results which uses a total of twenty (5N) function 
evaluations for a block of four steps (N+4) . Thus a sixth order 
Runge-Kutta is gotten which uses live function evaluations per step. 
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The above discussion illustrates the flexibility of these "Block" 
Runge-Kutta's — high order of accuracy for the number of function evalua- 
tions, and variable order of accuracy possible. 

Rosser 63 has derived several such methods where the number of 
function evaluations are further reduced by curtailing some of the 
early steps, and accelerat&hg the convergence of the later steps. In 
these methods order of accuracy is changeable, and is changed by 
varying N or r. An error estimate is also given, and a procedure for 
changing N and h from block to block outlined. In general if N = 2r 
or 2r + 1 a Runge-Kutta of order 2r + 1 results. 

The Block approach improves conventional Runge-Kutta's and 
makes them comparable to predictor-correctors on speed. But the 
"Block" approach is applicable to predictor-correctors also and will 
improve them too. 

Rosser has derived methods up to N - 8, i.e., up to a tenth order 

(a) 

Runge-Kutta. In the following formulae y is used in the error 
,s 

estimates and is — ■ evaluated in [x^, x^]. 
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One-point formulas: 


2ft - yo -hyc + 


h*y w 



Two-point formulas: 

yi - i/o - g (2/0 + yi ) , 

jj <») 

V* ~ I/o = 2Ay x ' + 2|-. 

Three-point formulas: 

i <. < , 4 » , \i A 4 y (4> 

2/1 - - p ( a 2/o + Slh — J/s ) H gp , 

I/s J/o = g (i/o + 4jn 4* j/j ) gQ- , 

2 /» — j/o = (yo + 3 j/j ) + — g — . 
Four-point formulas: 

I f» “ 3f» “ 55 - 5y s ' + y/) - 

I/s - yo - ^ (yo' + 4y»' + y*') - , 

ys — yo 53 ^ (yo' + 3y/ + 3 y/ + y*') — , 


y« - y« = y (2yx' - y/ + 2y s ') + 




y« - y* = p dOyo' - Sly?' 4 136y,' + 31y,') - 


351AV" 
160 * 


Five-point formulas: 

yi - yo = — (251y 0 ' + 646y»' - 264y,' + 106y,' - 19y,') + 
y* — yo = (29yc' + 124y x ' + 24y/ + 4y*' — y/) + , 

yj — yo = jsQ (Dy 0 + 34y/ + 24y/ + 14y»' — jft ) + , 

y* - yo = ~ (7yo' + 32y/ + 12y»' + 32y s ' + 7y/) - , 


y* - Vo - ~ (1<W - ifV + 120ys' - 70y»' + S5y,'> + , 

y* - y» --- ”j ( j(32t)y,' + l.Vlily/ - 12 1 Gy,' + 4'»9y/ - 31y.') + 


A. 126 


A. 128 


A. 129 


A. 130 
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* “ W “ 555 + 45 V - 1216 y/ + 1590 m' + 329y/) 

_ 1«V W 

45 ’ 

^ ~ 2/0 = 4) (lW + 8ly/ — 64y/ 4- Sly/ + lly/) — . 

Six-point formulas: 

* “ * “ 1^0 ^l' 9 * + 1427 ^ “ 79S ^ + 482y/ - 173y/ + 27y/) 

_ 863A 7 y 0) 

€0480 ’ 

V* " Vt = 4 (28y#r + my ' + 14y * + 14 ^' “ 6y/ + y/) 

_ 37Ay 7> 

3780 ’ 

» "* » 38 ^5 (17m + 73 y/ + 38y/ + 38y»' - 7y/ + y/) 

my y 

2240 ’ 

2A 

y« — y« = ^ (7y/ + 32y/ + 12y/ + 32y/ + 7y/) 

_ 8A 7 y 0 ‘ 

945 ’ 

V* ~ y ° * J| ( l9 »' + '5y/ + 30y/ + 50y/ + 75 y/ + i9y/) 

_ 275A 7 y <7) 

12096 ' 

ot 

M — M = jq 01 V\ ~ 14y/ + 26y/ - 14y/ + 11 y/) 

. 4lA 7 y 0) 

140 ‘ 

Seven-point formulas: 

W ~ “ ^so( 19087 m' + 651 12y/ - 4G461y/ + 37504y/ 

- 2021 ly/ + G3l2y/ - S63y/) + ~iC , 

24192 

V* - V* ■ (1139y/ + 5640y,' + 33 m' + 1328y/ 

- S07y/ + 264y/ - 37 m' ) + , 


A. 131 


A. 132 
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* — * — ( 0S5* + 3240* + 1161*' + 2176*' 


- 729 *' + 216*' - 29*') + , 

Mb 


y* — y o = ( 14'V + 696*' + 192*' + 752*' 


+ S7*' + 24*' - 4*') + 


y> ~ yo= j7~. (743* + 34SO*' + 1275*' + 3200*' 


SftV” 

945 


L S (8) 


+ 2325*' + 112$*' - 55*') + — - V , 

y 24192 ’ 


y* - u*- (41*' + 216*' + 27*' + 272*' 


+ 27*' + 216*' +41*')-^*-, 

1400 


y* V» — ^40 (731* — 840* + 8547*' — 1164S*' 


+ 14637*' - 7224*' + 4417*') + 


y *~ Vo= ( - myo> + 1312yj> “ 2048*' + 3132*' 


Eight-point formulas: 


- 2048*' + 1312*' + 115*') - mtA>y< * > 

14175 * 


y* - » - 120960 (30799 ^' + 1 39349*' - 121797*' + 123133* A. 133 

- 8S547*' + 41499*' - 11351*' + 1375 *' ) - 

362SS00 ’ 

* - * = 3^0 (1107*' + 5SG4*' - 639*' + 244$*' 


- 1927*' + 936*' - 261*' + 32* 1 - 

16200 * 


V*~ U»= (1325* + 6795*' + 1377*' + 5927* 


- 3033*' + 1377*' - 373*' + 45* > - 


o m*y w 
44800 ’ 


2 A 

y*- V 0 = g-jT ( 139*' + 724*' + 108*' + S92j /,' 


- 53*' + 108*' - 32*' + 4* > - 

14175 ’ 
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2/3 - 2/o = j (S341j/ 0 ' 4 46030y/ + lolOy/ 4 63670y/ 

— 800y/ + 341S6y/ — 9830 y/ 4 229QyJ — 245 y%) 

25 A V" 

+ 3584 ’ 

y* - yo= (401y 0 ' 4- 2232?// 4 18 y/ 4 3224y,' 

— 360y 4 4 2664i/j 4 lSSy#* 4 72y/ — * 9y/) 

9 h l0 y m 
+ 1400 ' 

- y°= (21361y 0 / 4 116662y/ 4 6958y/ 4 155134y/ 

4 7S40y 4 ' 4 105154y/ 4 7457Sy/ 4 31882^ — 1169y/) 

8183AV 10) 

+ 1036800 ' 

y* - Vo = (989yo 4 58 i$Sy/ - 92Sy/ 4 10496y/ 

~ 4540y/ 4 10496y s ' - 92Sy/ 4 5S8Sy/ 4 9S9y/) 

_ 2368&V 10 

467775 * 
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:>h 


yi — *o ™ ^ 1431*o' + 7345*/ + 1395*/ + S325*/ 

+ 27252// + 3411i// ~ 495*/ + 55*/) - 


Zfc — * yo — j^q (41t/o + 216y/ + 27y/ + 272ya # 


y? - y» - j^Q (751*/ + 3577*/ + 1323*/ + 2989*/ 


+ 27*/ + 216*/ + 41*/) - 


+ 2989*/ + 1323*/ + 3577*/ + 751*/) - , 


Ql. 

*» — *«= ^ (460*/ — 954*/ + 2196*/ — 2459*/ 

+ 2196*/ - 954*/ + 460*/) + 


3956AV” 
14175 ‘ 


Xine-pomt formulas: 

h 


*i-*»= 3g28i0O (107001 'Sfo + 4467094*/ - 4604594*,' 

+ 5595358*/ - 5033120*/ + 314633S*/ - 1291214*/ 


+ SWI»'-^') + «gg£. 


*i - *0 


113400 


(32377*/ + 182584y/ - 42494*,' 


+ 1200SS*/ - 116120*/ + 74728*/ - 31154*/ 

n .tn (io> 

+ 7624*/ - 833*,') + 


Vi-yo= (12881*/ + 70902*/ + 343S*/ + 79934*/ 

- 56160*/ + 34434*/ - 14062*/ + 3402*/ - 369*/) 

25 AV 10 ’ 

**■ 35S4 ’ 


V* - y« - Jfifi (4063*/ + 22576*/ + 244*/ + 32752*/ 

- 9080*/ + 9232*/ - 3956*/ + 976*/ - 107*/) 

, 94A i y , “ > 

+ 14175 ’ 


A. 134 
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Butcher [56,57] also creates hybrid methods which are a combination 
of RKMS's and predictor-corrector methods. Butcher [57] lists RKMS's 
with = 2,3, and 4 (see section A. 4) and of order p = 2S n + 2, with 
V — 4 • Methods up to = 15 are proved to exist. On basis of various 
tests these methods are highly recommended (they compare favorably 
with standard predictor-correctors). In the opinion of the author, 
the Rosser "Block" methods are, on the whole, better than these. 

Gear [64], Gragg arid Stetter [65], and Dahlquist [66] present 
methods similar to Butcher's methods. 

Another special class of Runge-Kuttas is given in Gates [26]. 

Gates formulates explicit Runge-Kuttas in which each stage (each k^) 
is independent of the other stages. In the conventional Runge-Kuttas 
each stage is not independent of the other stages. Each independent 
stage of the Gates formulae corresponds to a standard quadrative 
formula, e.g. , Gauss, Radau, Lobatto, etc. Making each stage inde- 
pendent of the others has a penalty — V increases for a given order 
p as compared to the usual Runge-Kuttas. For example, for p = 4,5, 
and 6 Gates forms have V = 5,7 and 11 respectively. The advantage 
of Gates' formulation is the flexibility gotten by having the co- 
efficients of one stage independent of another; another advantage is 
that deriving equations for these methods are much simpler than the 
usual RKE's and. so high order methods can be derived fairly easily. 

Stoller and Morrison [67] and Day [68] discuss similar quadrative 
Runge-Kuttas. 

From the computational time point of view these methods are ob- 


viously not recommended. 
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f . 6 General Comments for Further Development of the Runge-Kutta 

Class of Methods : 

The previous sections have tried to sample various types of Runge- 
Kuttas available in current literature. 

The field of Runge-Kuttas is wide open and improvements are 
possible. Runge-Kuttas in general arrive at an accurate solution at 
the end of an integration step, by combining solutions (k^’s) at 
intermediate steps which are not, in themselves so accurate. So one 
way of improving Runge-Kuttas would be to choose the coefficients 
of each stage in such a way that the k^'s would be accurate estimates 
of f (T,Y) . 

This should be possible without significantly increasing V, 
unlike the Gates [26] approach (see section A. 5). Then as each k i 
would be an accurate estimate of f (T,Y) a hermitian curve fit would 
be possible giving a higher order of accuracy. 

The most obvious improvement needed by Runge-Kuttas is to reduce 
their stages V per step, this would make them comparable to predictor- 
correctors on computational time. In all other aspects (except error 
estimation), i.e. , stability and self-starting features Runge-Kuttas 
are superior to predictor-correctors. So the second area of improve- 
ment lies in reducing V. This is possible to some extent by using 

i 

the "Back Step" and "Front Step" methods outlined in section A. 1-2, 
or by using the Shanks' technique outlined in section A. 1-5. Con- 
ventional RKE's always have c^=0 and usually have a c^=l , hence the 
corresponding to the c^l of the (T ,Y n ) step could be used for 
the k^ of the (T^ step reducing V by one. This would cost 

something — accuracy, because the solution Y , available at the end of 
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(T ,Y ) gives a more accurate estimate for k. than is possible by using 
the k^ of the previous step. Also if two or more c/s have the same 
value, then only one k^ corresponding to the first of these c^s need 
be evaluated. For the rest of these k^*s no more function evaluations 
are necessary, they can use the same value as the first of these k^'s. 
Again V is reduced but a loss in accuracy would occur. For example, 
a RKE(4,4) which has c^=0, Cj® y, c 3 = \* c 4 *1 would become a RKE(4,2) 
by using the two techniques outlined above — k^ at (I ^ t Y ^ would 
be the k y of (T ,Y ), k_ at (T . - ,Y ,,) would be k„ of the same step. 

It is not possible to say, off hand, whether the "Back Step" or "Front 
Step" technique would be more accurate than this one of interchanging 
k/s. Obviously the interchange techniques could also be used to 
further reduce the V of a "Back Step" or "Front Step" method, with 
(probably) a further loss of accuracy, or the Rosser [63] (see section 
A. 5) "Block" technique could be used to further speed up these methods. 

A third area of development (not speed improvement) lies in the 

f 

RKSI class of methods. Here a calculation procedure for — or a 
Jacobian is required. As this is available it might as well be used 
to improve the order of accuracy, though at the cost of extra compu- 
tational effort. 


n+1 


V 

l 


Y + i w.k, + ; u.g. 
n i£i 1 1 i=i 11 


V 

l 


k. = hf (T + c h, Y + I a.,k.) 
i n i n ^ lj y 

H 

‘ h2f *< T n + d i h > Y n + £ Vj> 

L. < V 

l 


A. 135 
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The Rosser [63] (see section A. 5) should be applicable to speed up 
any type of Runge-Kutta. 

Many authors have created hybrid methods which combine Runge-Kutta 's 
with predictor-correctors, to try and combine the advantages of both 
methods (see section A. 5). Similarly hybrid methods combining Runge- 
Kuttas with extrapolation and other methods should be possible. 

The Rosser [63] "Block" technique (see section A. 5) could be 
adapted to further speed up these hybrid methods. These hybrid 
methods may grow to such forms that Runge or Kutta would hardly 
recognize their method. 

There are other areas in which Runge-Kutta' s can be developed 
and improved, it is hoped that present and future mathematicians will 
do so. 
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APPENDIX II 
MULTI-STEP METHODS 


This appendix lists many of the multi-step methods that are in the 
current literature. 

The most common form of multi-step methods is 


I <«.y 4.1 . + h &. y' .) = 0 . 
, n 1 n+l-i 1 ■'n+l-i 

i=0 


When Sq = 0, then the equation becomes an explicit equation for y n+ ^; 
such equations are called Predictor (P) equations. When 6^ ^ 0, 
the equation becomes an implicit equation for y n+ ^j i.e., the equation 
must be iterated for y n+ ^* Such equations are called Corrector (C) 
equations. 

Although the Predictor (P) equation can be used with one derivative 
evaluation (E) at each step of the numerical solution, usually a combina- 
tion of the P and C equations is used to solve the initial value prob- 
lem. The various modes of these P-C combinations are PECE, PEC, PE(CE) S 
s 

and P(EC) . The first two modes are not iterative methods. First, 

y is predicted. Then, there is a derivative evaluation of y' „ . 
n-i-± n+ 1 

Then y n+1 is corrected using y^ + ^. For the first mode a new y^ is 
evaluated. The third and fourth modes simply iterate s times on the 
corrector. 

The Predictor and Corrector equations presented below are arranged 
in related groups. The format of each equation is in the following 
order: 1) the equation number S, 2) the order m and number of 

backsteps n, written (m,n), and 3) the a and 3 coefficients. The 
general format is 
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S. (m,n) (—-) (a. + a, + • •• + a ) 

a o 1 2 n 

+ d=-) (8 n + B. + ♦•• + B ). 

B 0 0 1 n 

Euler Predictor: [556] 

1. (1,D (j) CD + (~) (0 + 1) 

Milne Predictor: [556] 

2. (4,3) (j)(0 + 0 + 0 + 1) + (0 + 8 - 4 + 8) 

Mi llman-Kl op f enstein Predictor : [277] 

3. (4,4) (p> (-.29 - 15.39 + 12.13 + 4.55) + (y) (0 + 2.27 + 6.65 

+ 13.91 + 0.69) 

Craine-Klopf enstein Predictor: [95j 

4. (4,4) (y) (1.547 - 1.867 + 2.017 - 0.6973) + (y)(0 + 2.002 - 2.031 

+ 1.818 - 0.7143) 

Hermite Extrapolation Predictors: [556] 


5. 

(3,2) 

(y) (-4 + 5) + (y) (0 + 4 + 2) 

6. 

(4,3) 

(y)(-9 + 9 + 1) + (y) (0 + 6 + 6) 

7. 

(5,3) 

(y) (-18 + 9 + 10) + (y)(0 + 9 + 18 + 3) 

8. 

(7,4) 

(■—) (-128 - 108 + 198 + 47) + (y) (0 + 16 + 72 + 48 + 4) 

9. 

(9,5) 

(|) (-475 - 1400 + 600 + 1150 + 131) + (y) (0 + 25 + 200 


+ 300 + 100 + 5) 
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Nystrom Predictors: [556] 

10. (2,1) (y)(0 + 1) + (y)(0 + 2) 

11. (3,2) (y)(0 + 1) + (y)(0 + 2 + 0) 

12. (4,3) (y)(0 + 1) + (y)<0 + 7 - 2 + 1) 

13. (5,4) (y)(0 + 1) + (|)(0 + 8 - 5 + 4 - 1) 

14. (6,5) (y)(0 + 1) + C^)(0 + 269 - 266 + 294 - 146 + 29) 

Same Predictors from method of undetermined coefficient: [556] 

15. (5,4) (y)(0 + 0 + 1) + |(0 + 21 - 9 + 15 - 3) 

16. (5,4) (y) (0 + 0 + 0 + 1) + -^(0 + 9 - 32 + 64) 

17. (5,3) (y)(-9 + 9 + 1) + (y)(0 + 6 + 6) 

18. (5,3) (y)(-8 + 9) + (j)(0 +17+14-1) 

19. (5,3) (y) (-7 + 9 - 1) + (y)(0 + 16 + 10 - 2) 

20. (4,3) (y) (-54 + 45 + 10) + (y) (0 + 24 + 42) 

21. (4,2) (y) (-4 + 5) + (y)(0 + 4 + 2) 

22. (5,4) (y) (—1 + 0 + 1+1) + (y) (0 + 3 + 0 + 3) 

Schoen Predictors: [461] 

23. (5,5) (y) (-0.01745885 + 1.29864818 - 0.13318934 - .14799999) 

+ (y) (0 + 3.0001344 - 3.05980708 + 2.98844518 - 1.66187410 


+ 0.32137111) 
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24. (6,6) (y)(-. 92790378 + 2.05094161 + .54664912 - 0.67100875 

+ 0.0005296 + 0.00126884) + (y) (0 + 3.60328436 

- 3.60284254 + 5.40345116 - 4.77344249 + 1.80481406 

- 0.29749495) 

Krough Predictors: [556] 

25. (4,4) (A) (y)(l + 1) + (y^)(0 + 119 - 99 + 69 - 17) 

26. (4,4) (B) (y)(4 + 3) + (yy)(0 + 103 - 88 + 61 - 15) 

27. (5,5) (yy)(-l + 32) + (yy-y)(0 + 22,321 - 21,774 + 24,216 

- 12,034 + 2391) 

28. (6,6) (_L) (-11 + 23) + ( 17) 1 28Q )(0 + 62,248 - 62,255 + 101,430 

- 76,490 + 30,545 - 5079) 

29. (7,7) (jt)(-21 + 31 > + < 6047800 ) ( ° + 2 ’ 578 ’ 907 ~ 2 «^4,408 

+ 5,615,199 - 5,719,936 + 3,444,849 - 1,149,048 
+ 164,117) 

Adams-Bashforth Predictors: [556] 

30. (2,2) (y) (1) + (y) (0 + 3 - 1) 

31. (3,3) (y) (1) + (-^)(0 +23-16 + 5) 

32. (4,4) (y) (1) + (y^)(0 + 55 - 59 + 37 - 9) 

33. (5,5) (y) (1) + (y|o )(0 + 1901 " 2774 + 2616 " 1274 + 251 > 

34. (6,6) (y) (1) + (i^q)(0 + 4277 - 7923 + 9982 - 7298 + 2887 - 475) 
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35. (7,7) (j) (1) + ( 6Q ^ 8Q ) (0 + 198,721 - 447,288 + 705,549 - 688,256 

+ 407,139 - 134,472 + 19,087) 

36. ( 8 , 8 ) (j) (1) + (j20796O )(O + 434 > 241 ~ 1.162,169 + 2,183,877 

- 2,664,477 + 2,102,243 - 1,041,723 + 295,767 

- 36,799) 

Milne Corrector: [556] 

37. (4,2) (i )(0 + 1) + (|)(1 + 4 + 1) 

Hamming Corrector: [556] 

38. (4,2) (|)(9 + 0 - 1) + (|)(3 + 6 - 3) 

Milne-Reynolds Corrector: [361] 

39. (5,3) (-|)(1 + 7) + (j|2 )(65 + 243 + 61 + 1) 

Correctors from method of undetermined coefficients: [556] 

40. (4,2) (i)(4 + 1) + (|)(2 + 4) 

41. (5,3) (y)(0 + 0 + 1) + (|)(3 + 9 + 9 + 3) 

42. (5,3) (|)(1 t 1 + 1) + ( 3 ^) (13 + 39 + 15 + 5) 

43. (5,3) (-|)(1 + 1) + (~r)(17 + 51 + 3 + 1 ) 

44. (5,3) (-|)(2 + 1) + (y^) (25 + 91 + 43 + 9) 

45. (5,3) (-~)(9 + 9 - 1) + <jj )(6 + 18) 

46. (5,3) (|)(9 + 1 - 1 ) + (^0(10 + 22 - 8 ) 
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47. 

(5,3) 

(y) (9 - 1 - 1) + (yy) (8 + 14 - 10) 


48. 

(5,3) 

(yy)(45 - 9 - 5) + (yy) (12 + 18 + 18) 


49. 

(5,3) 

(~)(9 - 3 - 1) + (y)(2 + 2 - 4) 


Norsett A(a) 

stable Correctors: [389] 


50. 

(4,3) 

(yy)(18 - 9 + 2) + (yy)(6) 


51. 

(5,4) 

(yy)(48 - 36 + 16 - 3) + (yy)(12) 


52. 

(6,5) 

(y|y)(300 - 300 + 200 - 75 + 12) + (yyy) (60) 


53. 

(7,6) 

(— i 2. 

v 147)(360 - 450 + 400 - 225 + 72 - 10) + (— 

:) (60) 

Adams 

-Moulton Correctors: [422] 


54. ' 

(2,1) 

(y)(l> + (f)d + 1) 


55. 

(3,2) 

(y) (!) + (yy) (5 + 8 - 1) 


56. 

(4,3) 

(y) (1) + (yy)(9 + 19 - 5 + 1) 


57. 

(5,4) 

(y) (i) + (y|y)(251 + 646 - 264 + 106 - 19) 


58. 

(6,5) 

(y) (D + (y^o) (^75 + 1427 - 798 + 482 - 173 

+ 27) 

59. 

(7,6) 

(y) (!) + ( 6 o; 48 o )(19 » 087 + 65 ’ 112 ” *6,461 + 

37,504 


- 20,211 + 6312 - 863) 

60. (8,7) (y)(l) + ( 1207960 ) (36*799 + 139,849 - 121,797 + 123,133 

- 88,547 + 41,499 - 11,351 + 1375) 
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61. (9,8) (j) (1) + ( 3 >6 28 > 800 )(1>070>017 + 4 > 467 > 094 " 4 , 604, 594 

+ 5,593,358 - 5,033,120 + 3,146,338 - 1,291,214 
+ 312,874 - 339,533) 

62. (10,9) (y)(l) + ( 7> 257 > 600 )(2?082>753 + 2 > 449 » 717 “ 11,271,304 

+ 16,002,320 - 17,283,646 + 13,510,082 - 7,394,032 
+ 2,687,864 - 583,435 + 57,281) 


Fehlberg Correctors: (422] 


63. 

(3,2) 

(y)(4 + 1) + (j)(2 + 4) 

64. 

(4,3) 

(jj) (9 + 9 - 1) + (jj) (6 + 18) 

65. 

(5,4) 

(yy)(16 + 0 + 0 + 11) + (~)(3 + 10 + 0 + 6 + 1) 

66. 

(6,5) 

( 21, 1 319 ) ( ° + 0 + 13,300 + 8775 - 756) + ( 21 3 19 ) (672 ° 
+29,700 + 13,500 + 21,300) 

67. 

(7,6) 

( lU60 )(243 + 0+ 125 + 0 + 0 + 632) + (^~) (120 + 567 



+ 0 + 600 + 0 + 405 + 72) 

68. 

(8,7) 

( 82,49 1 Q,048 )(0 + 0 + 28 »1° 7 > 625 + 0 + 29,545,048 


+ 2,783,200 - 1,994,625) 

+ ( ~g2 ' t 4$o , Q43 } C 2 4 , 2 9 9 , 5 20 + 126,015,750 
+ 21,168,000 + 0 + 70,634,970) 


Wesson Correctors: [422] 

69. (3,2) (-|)(1 + 1) + (|)(3 + 8 + 1) 

70. (3,2) (~)(9 + 16) + (^>(109 + 328 + 55) 
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71. (5,4) 

72. (5,4) 

73. (5,5) 

74. ( 8 , 8 ) 


Rodabaugh-Wesson Correctors: [422] 

75. (5,4) <-£)(l + 2 + 4 + 9 ) + ^q ) (3703 + 15,518 + 6168 

+ 10,898 + 1873) 

76. (7,6) (-^-Hl + 2 + 4 + 8 + 16 + 33) + (^ 3 —^) (128, 627 

+ 642,168 + 130,167 + 693,632 + 143,137 + 399,240 
+ 61,469) 

77. ( 8 , 8 ) (i)(l) + ( 37 6 - 2 | ;8Q 0> d, HI >267 + 413,709 - 3,449,594 

+ 3,285,358 - 2,145,620 + 836,338 - 136,214 
- 17,126 + 7297) 

78. (9,8) ( 2 56 q -- q I 6 )(9784 + 20,133 + 41,040 + 79,775 + 159,816 

+ 319,691 + 639,792 + 1,289,985) 

+ ( 2,56i,016 )(725 - 340 + <l - 150 - 740 
- 280,710 + 6,541,620 - 1,808,250 

+ 5,630,940 + 244,290 + 2,458,620 


(y^)(l + 0 + 0 + 15) + ( I f J 2 q) (3611 + 16,006 + 5496 
+ 15,466 + 3341) 

(-jj )(8 + 2 + 1 + 5) + (^ 4 ) (767 + 2638 + 168 + 1282 + 185) 

( I )(1) + ( 2,440,080 )(859 ’ 838 + 2 ’ 143 ’ 299 “ 80 2 >706 

+ 267,244 - 18,396 - 9199) 

(y)(l) + ( 3 ~ ,628,80d )(1 ’ 147)591 + 3 ’ 846 ’ 502 " 2,432,522 
+ 1,251,214 + 397,060 - 1,197,806 + 880,858 

- 307,718 + 43,621) 


+ 345,330) 



